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In recent years regional fine mesh models have gained an important 
place alongwith global models both in operational numerical weather predic- 
tion and meteorological research. Compared to global models, the horizontal 
resolution in these models can be increased with ease, to provide detailed 
numerical guidance for the region of interest. Also these models are parti- 
cularly relevant to work on the development of physical parameterization 
schemes. 

The motivation for this study is the development of one such regional 
weather prediction model with an emphasis on certain aspects of the land 
surface processes. The thesis has seven chapters. Salient features of various 
regional models currently in use all over the world for operational and rese- 
arch purposes are briefly reviewed in the introductory chapter. This provides 
a perspective to the performance of the model discussed in this study. 
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The second chapter contains a detailed description of the dynamical frame- 
work of the model. The domain chosen for the model is in the summer 
monsoon region from 45°E to 114.4°E West-east and 9.4°S to 39.4°N South- 
north. The horizontal resolution for control runs is 1° lat/long which is 
reduced to 0.5° for sensitivity experiments. In the vertical the model uses 
the sigma (cr = p/p^ ; p g E surface pressure) coordinate system in which 
the earth's surface with orographic features is a coordinate surface 0=1. 
The non-linear horizontal advection is computed by the semi-La grangiam 
advection scheme that controls non-linear computational instability. The 

time-differencing is based on a semi-implicit breakdown in which the linear 
terms responsible for the excitation of gravity waves only are treated impli- 
citly while the slower Rossby modes are handled explicitly. The semi- 
implicit casting of the time-differencing is known to provide computational 
efficiency over a fully explicit treatment. An Arakawa C grid for the 

horizontal lattice structure ensures rapid adjustment between mass and 

flow fields i.e. geostrophic balance in higher latitudes. The forecast variables 

are not determined directly but through a model defined variable whose 
evaluation requires the solution of a three dimensional Helmholtz equation. 
Using a similarity transformation this three dimensional equation is reduced 
to a system of two dimensional equations which is solved by the fast fourier 
transform method. The model uses a dynamic initialization scheme. This 
scheme involves one time step forward integration of the model with all 
the physical processes included. It is followed by a number of cycles of 
forward-backward integration while keeping the non-linear tendencies fixed. 
Two options of fixed and time dependent lateral boundary conditions are 
currently available in the model. 
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The physical processes in the model are described in the third chapter. 
They include a comprehensive state of the art parameterization of the 
surface fluxes through the similarity theory, definition of a planetary boun- 
dary layer, cumulus convection through a version of Kuo's method, a band 
model radiative transfer and cloud model radiative feedback processes, 
solution of surface energy balance equation and the use of the envelope 
orography. 

Details of the case study chosen for the control and sensitivity experi- 
ments are given in chapter 4. The landfall of a tropical cyclone containing 
spiral rainbands in its cloud cover and precipitation patterns is chosen for 
the study. This cyclone formed in the southern Bay of Bengal during early 
May 1979, moved northwest to the east coast of southern India where it 
made its landfall around 12 May 1979. The model was initialized with the 
FGGE III-b data sets for 11 May 1979 00 UTC, and integrated for 72 
hours. The forecast circulation features compared very well with the obser- 
vations. The storm track was realistic. However the rainfall rates were 
poorly predicted. 

One of the crucial parameters of the land surface processes in estimat- 
ing the surface fluxes of moisture in the model is the ground wetness/ surface 
soil moisture variable, which is empirically specified. This is recognized 
as an area of weakness in the current modelling effort, leading to erroneous 
forecast of rainfall rates. Chapter 5 highlights this problem and gives details 
of a new parameterization scheme based on carefully designed moisture 
budget studies. A detailed analysis of the new scheme is made and the 
results of its impact on the estimation of various surface parameters are 
compared with the old scheme. A major improvement in the surface moisture 
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flux estimation is noticed. The forecasts show improvement in the rainfall 
rates and spatial positions of rainfall regions compared to observations. 

However there is scope for improvement around the storm region. In the 
next chapter sensitivity experiments for evaluating the role of horizontal 
resolution are discussed. Two different horizontal resolutions ( l°xl° 
; 0.5°x0.5° ) are chosen for prediction with the old and new parameterization 
schemes for the ground wetness. These experiments are found useful in 
assessing the model capability in resolving the details of the complex geo- 
metry of the observed precipitation patterns from the new improved para- 
meterization scheme for the ground wetness. Considerable improvement 
was noticed in the forecasts when the resolution was halved to 0.5° in both 
the spatial positioning of the rainfall region and the amount of rainfall 
around the storm region. The final chapter summarises the important conclu- 
sions. 

Thus, the study shows the development of a regional numerical weather 
prediction model for use in the tropics. In a case study of the landfall of 
a tropical cyclone, the model predicted reasonably well the storm track 
and the rainbands of the cyclone as it made the landfall. With the incorpora- 
tion of a new parameterization scheme for the ground wetness, the model 
also predicted realistic rainfall rates over the storm region, at the same 
time enhancing the accuracy of the geometry of the spiral rainbands com- 
pared to observations. 



CHAPTER 1 


INTRODUCTION 


In the tropics mesoscale systems dominate the regional weather and 
are often part of larger systems like tropical storms. Hence mesoscale 
models are useful for tropical weather prediction. There are basically two 
types of mesoscale models. There are those that make use of horizontal 
resolutions of the order of 50 km or larger and can resolve mesoscale pheno- 
menon whose scale is of the order of a few hundred km, these are quasi- 
static models. There is also a considerable current research interest in 
nonhydrostatic mesoscale models that use a horizontal resolution of less 
than 5 km. The latter are usually designed to study the life cycle and detai- 
led mechanisms of individual mesoscale systems. The present study covers 
the former class of quasistatic regional mesoscale models. These models 
are usually as sophisticated as any of the well known global models in terms 
of their handling of dynamical and thermodynamical processes. The global 
models are evolving very rapidly with a substantial increase in horizontal 
resolution. The performance of the high resolution global model has reached 
a level where one is frequently asking whether a very high resolution global 
model might be the appropriate mesoscale (quasistatic type) model since 
it doesn't have to face the problem of lateral boundary conditions. In the 
coming years perhaps we would be able to answer this question. 


Several promising mesoscale regional models have emerged in 
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recent years that deal with the prediction of heavy rainfall. Table 1 

contains a brief summary of some of these regional models. This 

* 

table, in three parts provides an outline of the NCAR/Penn State, 
JMA, NMC, French, ECMWF regional models and the FSU regional model 
presented in this thesis. This summary includes the horizontal and 
vertical discretization, time differencing schemes, the lateral boun- 
dary conditions, initialization, the diffusion terms and physical pro- 
cesses (surface processes, planetary boundary layer, cumulus 
parameterization, radiative transfer) and domain of application. 

While superficially, the mesoscale prediction models may appear 
as mere refined versions of limited-area, primitive equation, numeri- 
cal weather prediction models for the synoptic and large scale, the 
forecast accuracies of these models are however to a greater extent 
constrained by limitations imposed by (a) numerical difficulties 
related to finite differencing and lateral boundary conditions, (b) 
adequacy of parameterization of subgrid scale physical processes, (c) 
uncertainties in the model's initial conditions and (d) inherent pre- 
dictability of atmospheric mesoscale circulations. In the tropics 
invariably one of the major limitations of mesoscale models ini- 
tialized through synoptic and subsynoptic scale data is that they 
show limited skill in simulating the spatial and temporal distribution 
of precipitation fields despite well predicted flow patterns. This 
could be attributed to slow spin-up of mesoscale circulations in the 
model initialized with synoptic and subsynoptic-scale data and inade- 
quate representation of surface latent heat fluxes. In the present 

* 

List of acronyms is given on page xvi . 
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Model 

Horizontal & 
Vertical 
Discritization 

Time-Differencing 

Scheme 

Lateral 

Boundary 

Conditions 

NCAR/ 

Penn 

State 

limited 

area 

model 

The model variables are 
staggered in the verti- 
cal and horizontal. In 
horizontal Arakawa B 
grid (Arakawa & Lamb, 
1977)is used. In verti- 
cal 16 a levels are used 
where the resolution is 
higher in PBL . For the 
horizontal grid Polar 
stereo graphic, Lambert 
Conformal or Mercator 
map projection is used. 

Explicit scheme 
developed by Brown and 
Campana ( 1978) is used. 

A frequency filter 
(Asselin, 1972) is 
applied to all 
prognostic variables. 

i 

j Five types of lateral 
boundary conditions can 
be used. 

1) Fixed 

2) Time-dependent 

3) Time-dependent 
and inflow/outflow 
dependent 

4) Sponge 

5) Relaxation 

! 

JMA 

Tokyo 

Regional 

Spectral 

model 

Spectral representation 
in horizontal using a 
modified double Fourier 
series. The series 
composes of the orth- 
ogonal double Fourier 
series with mirror 
boundary conditions 
and additional bases to 
specify fluxes across 
the boundary. In verti- 
cal, finite differencing 
is carried out on 16 
equispaced o levels. 

Semi-implicit scheme 
is used for integrating 
the momentum, thermo- 
dynamic equation and 
continuity equation. 
Time integration of 
moisture field Is car- 
ried out by the leap 
frog scheme. 

Time-dependent b.c.s are 
prescribed by the global 
model forecast. Boundary 
relaxation is performed 
using a forward 
time-differencing scheme. 

NMC 

The 

Nested 

Grid 

Model 

In the forecast three 
grid resolutions are 
used. Stereographic 
coordinates and Eli- 
assen staggering is 
used in horizontal . 

In vertical, discriti- 
zation suggested by 
Arakawa (Arakawa et al . 
1974) is used. 

The time integration 
scheme is a two-step 
second-order Lax- 
Wendroff process 
(Phillips, 1962; Lax 
and Wendroff, 1964). It 
uses the Eliassen time 
staggering variables 
(Eliassen, 1956) . 

At the equator symmetric 
boundary conditions are 
chosen. 

' 

The 

French 

Weather 

Service 

Limited 

Area 

Model 

An Arakawa C-type grid 
! is used. Horizontal 
discritization is 
adapted from Sadourny 
(1975) and Burridge & 

• Haseler(1977) . In vert- 
ical 15 unevenly spaced 
a levels are used where 
resolution is higher 
in PBL. 

Leap-frog time integ- 
ration scheme is used. 
Linear part of the 
gravity term is treated 
semi-implicitly . 
Vertical advection and 
diffusion terms are 
treated implicitly. 

Fixed large-scale boundary 
conditions are used. 

Boundary relaxation 
proposed by Davies (1976) 
is used and is treated 
implicitly. 

ECMWF 

Limited 

Area 

Model 

| 

Fields are staggered in 
horizontal as well as in 
vertical. Arakawa c-type 
grid is used in hori- 
zontal . 

Semi-implicit time- 
stepping scheme is used 
following Robert et aJ. 
(1972). It is based on 
a time centered scheme. 

Time-dependent b.c.s 
derived either from a 
previous global forecast 
or from analyses given at 
intervals of a few hours. 
Boundary relaxation tech- 
nique suggested by Davies 
(1976) is used for smoothing 

FSU 

Limited 

Area 

Model 

Fields are staggered in 
horizontal as well as in 
vertical. Arakawa c-grid 
(Mesinger & Arakawa, 
1976) is used in hori- 
zontal . 

A one-step semi- 
Lagrangian advection 
(Mathur, 1983) and a 
semi-implicit time 
differencing scheme is 
used. 

Model has an option of two 
types of boundary conditions. 

1) Fixed boundary conditions 

2) Time-dependent 
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PHYSICAL 

PROCESSES 

Model 

Initialization 

Horizontal 

Diffusion 

Surface 

Parameters 

NCAR/ 

Penn 

State 

limited 

area 

model 

Vertical mode 
initialization scheme 
developed by Errico 
(1986) . 

A second order 
diffusion near lateral 
boundaries and fourth- 
order diffusion inside 
the domain is used. 

1 ■ 

Ground surface temperature is 
computed from a surface-energy 
budget following the "force 
restore" slab model developed 
by Blackadar (Zhang & Anthes, 
1982). Fixed ground moisture 
is used. Surface fluxes are 
calculated either Bulk- 
aerodynamical ly following 
Deardorff (1972) or using sim- 
larity theory. 

JMA 

Tokyo 

Regional 

Spectral 

model 

Nonlinear normal 
mode initialization 
is used. 

Horizontal diffusion is 
done by linear fourth 
order Laplacian. 

Daily analysed sea surface 
temperature. Ground surface 
temperature calculated with 

3 ground layers. Fixed ground 
moisture is used. Surface 
fluxes are calculated by simi- 
larity theory following the 
scheme given by Louis. 

NMC 

The 

Nested 

grid 

Model 

Initialized fields 
are extracted from 
globally nonlinear 
normal mode 
initialized data. 

None 

Ground surface temperature is 
not calculated. No ground 
hydrology is done. Heat and 
moisture surface fluxes are 
not calculated over land. Sur- 
face fluxes are calculated 
Bulk-aerodynamical ly. 

The 

French 

Weather 

Service 

Limited 

Area 

Model 

Non-linear normal 
mode initialization 
is used. 

Fourth order 
horizontal diffusion. 

ittT* ftmm 1 * tti 

i i » " TwiTStflflMt E ETtH# *^B' - J v-V ; 

; t 4t % 

ECMWF 

Limited 

Area 

Model 

Non-linear normal 
mode initialization. 

Linear fourth order 
horizontal diffusion. 

Ground surface temperature is 
computed by a slab model. 
Stability dependent surface 
fluxes are calculated using 
similarity theory (Louis, 

1979) . 

FSU 

Limited 

area 

Model 

Dynamic normal mode 
initialization 
proposed by Sugi 
(1986). 

Linear second order 
horizontal diffusion. 

Ground surface temperature is 
computed by surface energy 
balance . 

Surface fluxes are computed 
using similarity theory. 


1 In three cases namely stable, mechanically-driven turbulence and forced convection 
the prognostic variables above the surface layer are computed from k-theory. In the 
fourth case of free-convection above the surface layer a mixed layer is considered 
where the vertical mixing is visualized as taking place between the lowest layer and 
each layer in the mixed layer- Above this layer k-theory is used. 
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Model 

1 PHYSICAL PROCESSES 1 

| 

1 Domain 

PBL 

Cumulus 

Parameterization 

Radiation 

NCAR/ 

Penn 

State 

limited 

area 

model 

! 

Vertical eddy fluxes 
are either calculated 
Bulk-aerodynamical ly 
or using high reso- 
lution Blackadar PBL 
model (Zhang & Anthes, 
1982). The atmosphere 
stratification is 
divided into four 
categories based on 
bulk Richardson 
number . 1 

i 

A modified Kuo 
scheme (Anthes, 
1977) is used. 

Dry convective 
adjustment is 
also carried out. 

Short wave a long 
wave radiation par- 
ameterized. For 
short wave trans- 
missivity multiple 
reflection (Benja- 
min, 1983) is used. 
Clear-air transmis- 
sivities, back 
scattering coeffi- 
cients and longwave 
emissivity are 
functions of path 
length and precip- 
Itable water. Cloud 
effects are inclu- 
ded following Ben- 
jamin, 1983. 

Variable 

domain but poles 
| not included. 

I 

Spectral 

model 

ty and shear dependent 
eddy viscosity coeffi- 
cients are used. 

tion of rain is 
included. 

temperature. 


NMC 

The 

Nested 

grid 

Model 

Eddy fluxes in the 
vertical are computed 
using austausch form. 

A modified KUO 
scheme is used. 
Large scale sat- 
uration adjust- 
ment is done. 

Same as FSU. 

Hemispheric 
model . 

The 

French 

Weather 

Service 

Limited 

Area 

Model 

Vertical diffusion 
of fluxes is done 
using K-theory 
according to Louis 
et al. (1981) . 

Parameterization 
of deep convec- 
tion by Bougeault 
(1985) and Geleyn 
(1985) and super- 
saturation are 
included. Evap- 
oration of rain 
is also allowed. 

Simple radiative 
transfer scheme is 
used. Water vapor 
is the only active 
constituent for 
emission & absorp- 
tion. Emissivity 
and absorptivity 
are parameterized 
by Lepas et al . 
(1979). Effect of 
clouds is also 
included. 

Domain centered 
on Western 

Europe. 

ECMWF 

Limited 

Area 

Model 

Eddy fluxes are 
calculated on the 
basis of the mixing 
length theory. 
Stability and wind 
shear dependent 
diffusion 
coefficients are 
used. 

Kuo-74 cumulus 
model, supersat- 
uration large 
scale condensa- 
tion and non- 
convective cloud 
parameterization 
are included. 
Evaporation of 
rain is included. 

Radiation para- 
meterization is 
done following 
Geleyn and 
Hollingsworth 
(1979). Gaseous 
absorption of H 2 O 
COg and O 3 is in- 
cluded using band 
model theory. 

Cloudy sky is al- 
so parameterized. 

An arbitrary 
rectangle but 
cannot contain 
a pole or cover 
a complete zone. 

FSU 

Limited 

area 

Model 

For the vertical 
eddy fluxes there 
are two options: 

1 ) vertical 
structure based 
on GATE . 

2) K-theory, where 
diffusion 
coefficient 
depends on the 
mixing length, 
vertical wind 
shear and 
stability of the 
atmosphere. 

The model in- 
cludes parame- 
terization of 
deep moist con- 
vection following 
Krishnamurti et 
al. (1983a), 
shallow convec- 
tion (Tiedke, 
1984) and large- 
scale conden- 
sation. 

Long-wave 
radiation 
parameterization 
based on 

Harshvardhan and 
Corsetti (1984); 
includes effects 
of water vapor, 

C0 2 , O 3 and 
clouds. Short- 
wave radiation 
scheme Is similar 
to UCLA/GLAS GCM 
as described by 
Davies (1982). 

Variable 
domain 
but poles 
not 

included. 
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work we attempt to provide a realistic representation of the surface 
latent heat fluxes to overcome the aforementioned inadequacy through 
the development of a new ground wetness parameterization. 

The objectives of the present study are: 

(i) Development of a high resolution regional weather predic- 
tion model for the tropics. 

(ii) Assess the performance of the model in predicting the land- 
fall of a tropical storm. 

(iii) Formulate a procedure to find the initial values of ground 
wetness and thereafter develop a ground wetness parameteri- 
zation scheme which includes the effects of precipitation 
and other land-surface processes in its specification 
during the integration. 

(iv) Assess the impact of improved ground wetness parameteriza- 
tion in estimating the surface latent heat fluxes and sur- 
face temperature. 

(v) Conduct sensitivity experiments to see the impact of 
improved ground wetness parameterization on rainfall pre- 
diction for the same tropical storm as (ii) at resolutions 
(1° x 1°) and x V) . 

The present model is considerably enhanced in the dynamical, phy- 
sical and initialization aspects. Specifically this model has been 
designed for studies at very high resolution towards the examination 
of heavy rainfall events. The model utilizes the usual sigma frame as 
the vertical coordinate and permits the use of steep mountains. A one 
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step semi-Lagrangian advection (based on Mathur, 1983) and a 
semi-implicit time differencing scheme provide the basic framework for 
this model. The finite differencing grid utilizes the Arakawa c-grid 
(see Mesinger and Arakawa, 1976). This grid is staggered along the 
horizontal and the vertical coordinates. Fourth order finite dif- 
ferences are used in the formulation of most of the dynamics and phy- 
sics of the model. The domain chosen for the present study is -9.4°S 
to 39.4°N and 45°E to 114. 4°E. 

The tropical cyclone is an integral part of the weather activity 
in most tropical ocean basins. Its influence extends into the 
extratropics. Their evolution has been a subject of interest for 
decades. During the pre-monsoon period of 1979, when the cross 
equatorial flow was weak, a tropical depression in the Bay of Bengal 
on May 4, intensified into a hurricane on May 8, and eventually made 
landfall on the southeast coast of India causing a major disaster. 
Recently a numerical study of the genesis of this storm using a fine 
resolution (%° x J4°) regional model has been conducted by Tuleya 
(1988). It is of interest and importance to see the capability of the 
present regional model in forecasting the landfall of this tropical 
storm. 

An important area deserving further work relates to the proper 
treatment of fluxes of latent heat at the earth's surface. An impor- 
tant factor in the surface latent heat flux estimation is the parame- 
terization of ground wetness. In climate modelling significant impact 
of this parameter on the precipitation predictions has been 
demonstrated by many investigators such as: Yamakazi (1989), Rind 
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(1984) and several others. Nearly all the regional models outlined 
in table 1 including the present FSU regional model use very simple 
parameterization of the ground wetness. These formulations do not 
provide the changes in ground wetness due to precipitation and many 
other land surface processes. Therefore we propose a new ground wet- 
ness parameterization scheme to remove these deficiencies to some extent. 
This method makes use of observed measures of vertically integrated humidity 
budgets to provide a database for the surface humidity fluxes. These, in 
turn, are used to derive a statistical parameterization for the ground wetness. 
This parameterization is then used in very high resolution modelling experi- 
ments. 



CHAPTER 2 


DYNAMICAL FRAMEWORK OF THE REGIONAL MODEL 


2 . 1 Vertical and Horizontal Discretization: 


The dimensionless a coordinate is used in the vertical: 


a 



( 2 . 1 ) 


where p is the pressure and p s is the pressure at the earth's surface. 
A list of useful symbols is given on page xvii. 


The vertical extent of the model ranges from a=a^ at top of the 
model to ct= 1 at the earth's surface. There are N layers in the verti- 
cal. The thickness of each layer, measured by Act, is uniform. Fig. 1 
shows the vertical discretization of the model, in which, 


N = 9 and ct t = 0.1. (2.2) 

Variables u, v, q, and T are defined at the centre of a layers 
(dashed levels). The vertical velocity a and z are defined at the 
layer interfaces (solid levels). 


Arakawa C type grid is used for the horizontal discretization. 
This appears to be very well suited for the representation of 
divergence and for the formulation of a Helmholtz's equation which 
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Figure 1: Vertical discretization of the model. 
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arises in the use of the semi-implicit time integration scheme 
(McGregor and Leslie, 1977), 

The variable representation on the staggered grid are shown in 
fig. 2. Other thermodynamic variables such as q and T are defined on 
top of z. a and lnp s are also located at these points. 

2 . 2 Equations of motion : 

The governing equations for a semi -Lagrangi an formulation in the 
a coordinate system are written as follows: 


2a - -a fa ♦ Vlf^tan*) 
Dt oo a 


_3z _ RT 31nPs + Fi 


’3x 




r l 


g = -a |2 - u(f+— tan<f>) - - RT^ES + F r 

Dt 3a a 3y ^Ty V * 

Thermodynamic equation 

DT -aT 


kT fe + Dlnp s 
K [a "DT" 


Dt 3a 
Mass continuity equation 


H/ Cr 


Dlnpc 


3a 


- V*V 


Dt 3ct 

Moisture continuity equation 

M = + E 

Dt 3a 


(2.3) 

(2.4) 


(2.5) 


( 2 . 6 ) 


(2.7) 


Equations (2.3)-(2.7) form the prognostic equations of the model. 


To close the system of equations (2.3)-(2.7), the following 

diagnostic equations are used 

Hydrostatic equation 

RT Cp8 Bn 1 * 1 

ga g 3a 


3z 

3a 


( 2 . 8 ) 
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Figure 2: 


Horizontal discretization (Arakawa c type grid) . 
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Equation of state 


_ RT _ RT_ 

a ~ P ” <H?s (2.9) 

The quantity rr K in equation (2.8) is obtained from the Poisson 
equation : 

' K 


l-” K - 


1000 


( 2 . 10 ) 


The operator D_ in equations (2.3)-(2.7), which represents the hori- 
Dt 

zontal Lagrangian operator, can be expressed, in spherical polar coor- 
dinate system as 

D_0 _ 1_0 + u i_() v 3_() 

Dt 3t acos4> 9X a 34> (2.11) 

or, in Cartesian coordinate system as 


d_() _ a_() . u a_o 

Dt 3t 3x 


+ v 


i_() 

3y 


( 2 . 12 ) 


2.3 The Integration Scheme 

A one step, second order, semi -Lagrangian advection scheme is 
used to integrate the model. A special feature of this scheme is that 
it accounts for time changes of dependent variables and advecting 
velocity over the trajectory traced by the parcel in time step At. 
The inclusion of this formulation of the advective processes is 
expected to give better prediction of the phase and amplitude of 
rapidly developing disturbances as is generally the case in tropics. 
From a physical point of view it is a direct approach. To make this 
scheme economical this scheme is coupled to the semi-implicit time 
integration scheme. 


Equations (2.3) to (2.7) are transformed into a form suitable for 
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the semi-Lagrangian semi-implicit formalism using the following 
definitions : 


IfF cosi|> d\d(j) . 


(2.23) 


(2.14) 


II cos<J>dAd$ 

the horizontal average over the model domain. 

F' = F - F* 

the deviation of the variable from its mean. 

Combining the thermodynamic energy equation (2.5) and the mass con- 
tinuity equation (2.6), a model parameter P is defined as 

P = gz + RT* lnp s (2.15) 

so at the earth's surface, 


P s = gz s + RT lnp s 


(2.16) 


where z s and T s * are the surface topography and area averaged surface 
temperature, respectively. 

A static stability parameter y is defined by the relation 


Y _ E [»£X _ 9TI 
1 a [ a 8oJ 


a horizontal area averaged static stability is given by. 


* R 
^ a 


kT 

a 


8T 

3(7 


and the deviation is expressed by 

.* 


(2.17) 


(2.18) 


(2.19) 


y * = y - r’ 

Upon applying (2 . 13)- (2 . 19) to (2.3)-(2.5) and (2.7) the latter can be 
written as 

Ea , _i|a t V ( f,jtan4>) - gg - R(T*.T') HJPs + F X 

= + v(f + ~tan4>) - (gz+RT*lnp s ) - RT'ligs + F x 


or Du 3P 
Dt 3x 


+ v(f+-tan<J>) - RT 1 _^£g 
oCT a ox 


+ F x K Sj 


( 2 . 20 ) 



Similarly we obtain 


Dv + 9P 
Dt 6y 


DT 3T 
Dt + 0 3a 


-a— - 
da 

u(f+^tancj>) 

- kT 

£ + Dlnp s 

[a HF"J 


• 3T 1 
3a 


~w 


( 2 . 21 ) 


a + Dlnp s 
a DT" 


H 

~ Cn 


from a use of the mass continuity equation we obtain, 


D (T-tcT*lnp s ) 


v aa 

- + = o 


kT 1 _ 3T1 
a 3 a 


- kT' 

I 2 + v*v 


[3a J 


H 


D_ 

Dt 


[~4 


a D_ 
~R Dt 


g - KT*lnp s 
3P 


^ + cry lnp s - — = a 


y jro _ • [ kT ’ _ 3T' | 
R [ a 3aj 

kT' 3T 


KT' 


+ V*V + 


H 

<5 


* • 
y eg 


3a 


- kT' 


3a 
3a 

|2 + V-V + E- 
3a j 


or D 


1 3P 


Dt [y 3a 


alnp s + a = 0 


'a kRT’ 
ay^” 


3a 

3a 


+ V*V 


R H 


Cpcry 


s 3 


(2.22) 


Da = _3 + E _ 

TTk-t- -1 — 


Dt ~3a _ * (2.23) 

In equation (2.20)-(2.23) all the non-linear terms appear on the right 
hand side. 

A few more definitions are required for the semi-Lagrangian formula- 
tion, these are given below 

F = F(x,y,a,t), a function at any point x, y, a and at current 
time step t (2.24) 

F + = F(x,y,a, t+At) , a function at point x, y, a and at a future 


time step t+At 


(2.25) 


F° = F(x-a,y-b,a,t) , a function at time t at a point x-a, y-b 

from where the parcel arrives at x, y at 
time t+At (2.26) 

In semi-Lagrangian form DF/DT is expressed by 

DF d H F F + -F° 
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A semi -Lagrangian time average denoted by overbar is defined as 


p = F~ l ~ + F° 

2 (2.28) 
Upon application of the semi-Lagrangian - semi-implicit opera- 
tions the equations (2.20)-(2.23) and (2.6) are transformed into the 
following equations. 


^ + f = Sl 

dt 3x 1 
+ & = So 


dt 


3y 


djj [ 1 3P 

L r* 


+ alnp s 


I T 0 

, + <7 = So 


dt 


dH<J 

Tf 


= s 4 


(2.29) 

(2.30) 

(2.31) 

(2.32) 


and 


dH^Ps + ^ + M = o 
dt 3o 


(2.33) 

Here the forcing terms on the right hand side of equations (2.20)- 
(2.23) are replaced by their values at the origin. Equations (2.29)- 
(2.33) can be expressed as a set of 5 symbolic algebraic equations 






j = 1,2. 


(2.34) 


As described earlier, the superscript + refers to a future value at a 
grid point, o refers to the current value at the point of origin. The 
subscript j indicates equations corresponding to equations 
(2 . 29)-(2 . 33) . 

The individual terms in the equation (2.34) may be identified as 

At 3P 


Wj = u + 
*2 


v + 


2 3x 

At 3P 
2 3y 


(2.35) 


(2.36) 



17 


i ap 


y” da ”*“*'8 2 


w 3 
w 4 

w 5 = lnp s + 


+ CTlnp„ + 


At 


= q 


At 


v*v + 


do 

da 


Xh 


X-, = u 


At 

2 


3P 

ta " 2S 1 


x 2 

*3 


= V 


At 


8P 

2 [By S2 


7* Is * ‘"“Ps - (,j - 2s 3t 


x 4 

X5 = lnp s 


q + s 4 x At 

At 
2 


V*V + 


do 

da 


By eliminating various quantities from equations (2.35)-(2.39) 
equation for the model parameter P is obtained 

At 2 a fi* ap] At it aw 3 - Wc * 0 

T 7 p + W [7* 3crJ = 7 T v * w + ^ 5 U 


(2.37) 

(2.38) 

(2.39) 

(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 
a single 

(2.45) 


where W = i\ w i + 1^2 (2.46) 

Equation (2.45) is a three dimensional elliptic equation. This 
equation can be solved with appropriate boundary conditions for P. W + 
can be obtained from equation (2.34) making use of equations 
(2. 40) -(2. 44) . 

Equations (2.35)-(2.39) can then be used to give the prognostic vari- 
ables. At upper and lower boundaries the vertical velocity a=0 (see 
section 2.5). At upper boundary, therefore, (2.37) becomes 
1 3P 


y da 


= W 3 - cr T lnp s 


At lower boundary 
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1 3P 
y* 3 a 


= W 3 - alnp s 


where o=l 


From the vertical structure of the model, given in figure 1, we 
shall make the following definitions: 


NP1: Number of a layer boundaries 
cr^ (k=l (1 )NP1 ) : o at layer boundaries 
&k (k=l(l)N): a at layer means 


Acr k = CT k+l - a k k=l (1 )N 

(2.47) 

A9 k = ^k+1 “ 9 k k=l ( 1 )NMI 

(2.48) 

p k = e z k + RT k* ln Ps 

(2.49) 

P NP1 = e z s + RT s* ln Ps 

(2.50) 

lnp s = £NP1 _ g z s 

S RT S * RT S * 

(2.51) 


In writing the equations in finite difference form, the following 
boundary conditions are used (see section 2.5): 


Top of the model 

<7=<7t> <r=0 (2.52) 

Bottom of the model 

0=1 , 0=0 (2.53) 

The left hand side terms (Wj ) of equation (2.34) are written in finite 
difference form as follows: 

Starting from the expressions for and Wg from (2.35) and (2.36) 

For any layer k we can write, 

iEk+i + i£k 

3x 3x 


Wik = fi k + 


2 


(2.54) 
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and „ At 1&C + 1 + ^Zk 

w 2 k = v k + — ay ay 

L 2 j (2.55) 

(2.54) and (2.55) can be expressed as a single vector equation, i.e., 

— • ► —► <4 

w k = ix®lk + itfzk 


5 k * 4| vt V,i * 7p x 


(2.56) 


The divergence of (2.56) is expressed by, 

*4 . ,4 ♦ i§ [ j . 


(2.57) 


From (2.57) by replacing k by k-1 we obtain 

v 4 -i - * 4 -i * * [ ] 


(2.58) 


adding (2.57) and (2.58) and dividing by 2 we obtain, 

i ( v4 - v4-i> - + A| [ 1 2 V»Pk ^ Zlffcl ' 


For the top layer (k=l) equation (2.57) takes the form 


(2.59) 


V*Wi = V*Vi + 


- v#v l + 2 


At j V 2 P 2 + V 2 P X 


(2.60) 


For the lowest layer (k=N) equation (2.57) takes the form 


*4 - v4 ♦ 4| [ v * p NPi ; 7ap N ■ 


(2.61) 


We shall next examine the expression for W 3 : 


At general levels k and k-1 


1 [ p k+l " p k 


& k lnp s 


a k+l + 

. 2 “ 


w 3k-l 


.1 p k ~ p k 

\-l A(T k- 


— + frk-i ln Ps + — — 

1 L 2 


(2.62) 


(2.63) 


From equations (2.62) and (2.63) we obtain. 
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w 3k " w 3k-l 

jgr- 


L tk 


+ lnp s + 


p k+l " p kl 1 

* 


At 

2 


~^k . 

?k-l 

CT k+l " a k-l 


p k ~ p k-l 


Acr k-1 


]] 


TS57 


For the top layer equation (2.62) takes the form 


W 3 1 - 


*r 


p 2 “ P 1 


L<*2 “ CT 1 


+ & l lnp s + 


or W31 _ 1 




f p 2 " p l 1 

+ f a i + A ®i 

L ff 2 " °l J 

L “2“ J 

[ p 2 ~ p l' 

[ct 2 - aiJ 

+ ft- lnp s j 

Adi s + 1 


lnps + At °2 
2 




(2.64) 


(2.65) 


For the lowest layer we obtain, 


W3N = -1^ 1 — 

a %mi ?n a& nmi 


P NP1 “ P N 


ct NP 1 * ct N 
At general levels W5 is written as: 

f CT k+l ~ CT k 


„ lnp s _ I lnp * At a N 
«P1 2 inp s + 2 (2.66) 


»5k = ln Ps + ^ V-V k + ^ 


~^k 


£*, L . At v y At 

w 5k-l = lnp s + ~2 V * V k-l + “2 


£k~ g k-l 
Aa k-1 


(2.67) 

( 2 . 68 ) 


In the present case Ad k = Ad k _i , combining equations (2.67) and (2.68) 
we obtain, 


W 5 k + *5k-l _ .M (V-V k + V.V k _i) + At [<h< + i - <Hc-l | 

2 - tnp s + 2 2 2 ^ 2Ao£ J 

combining equations (2.59), (2.69) we obtain 


(2.69) 


*5k + »5k-i , lnPs * 41 (V.w k + 


At' 

2 

[V*P k+1 + 2V*P k + V*P k _j] 

2. 


4 J 


+ — [ g k+l ~ g k-l 

2 [ Sct£ 


(2.70) 
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For the top layer combining equations (2.60), (2.67) we obtain 

(2.71) 


W 5 l = ln Ps + 


[ v4 x - 


At 

2 


V*P 2 + yapj 
2 


+ At a 2 
2 Act 3 


or 


lnp s + 


— a 2 - w*-, - — 

2 SaJ W 51 2 


v*w n + 


At 

2 ' 

2 



V Z P 2 + 
2 


(2.72) 

For the lowest layer combining equations (2.61), (2.67) we obtain 

At 


%N ln Ps + 


V *% - *§ 


V^P NP3 + v 2 p n 
2 


2 So^ 


or 


lnp< 


— = ^5N _ At y.^, 


2 Act, 


N 


N 


At’ 

2 

vzp NPl + V * P N 

2 


2 J 


(2.73) 


(2.74) 


Upon elimination of various terms from the finite difference forms of 
the equations of motion, thermodynamic energy and mass continuity one 
obtains a single three-dimensional Helmholtz equation 


At 


yap + — 

V F 3 a [y" 3o 


1 3P 


■ ~z ™ - *3 ■ o 


(2.75) 

Its finite difference form is obtained by combining equations (2.64), 
and (2.70), i.e., 


i- [_y 


p k+l ~ p k 


AtT k 


k-1 


p k " p k 


^k- 


k-1 

1 J. 


2 f v * p k+l + 2VZp k + V2p k 


-1 


«3k " W 3 k-i 

_ [Wsk + %k-ll 


l j 

L 2 J 

(2.76) 


At 

2 

At 

2 


For the top layer, substituting for the last two terms in equation 

(2.65) from equation (2.72), we obtain, 

At v.w 3 W 31 _ W 51 

2 — zr soy ~z r 


« 1 _1_ f p 2 - p l 
?1 [a 2 - 


t ^ i,p s ♦ p§] 2 


(2.77) 


Substituting for lnp s in equation (2.77) from equation (2.51), we 
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obtain, 


1 _ 1 _ 
h* *»i 


p 2 


i a 2 ~ 


* oi P NP1 * [At] 2 fv zp 2 + V«Pil 

za^f rts» + [ 2 j i 4 j 


. At V-Wi W 33 _ W 51 gz s cr a 

2 2 T~ I&l (2.78) 

Similarly, at the lowest layer from equations (2.66), (2.74) and 
(2.51) we obtain 
1 1 


P NP1 

- P N ' 

i _ P NP1 

-1- + 

At 

2 

V ep NP i + V 2 P n 

L a NPl 

- ct nJ 

RT* 

Ad N 

2 


4 


^3N _ W 5 N At V*W N 

~S~ 2 2” 


gz S _L_ 

A % 

NP1 


(2.79) 


Thus equation (2.75) is written in a finite difference form as follows: 


1 

1 

p k+l p k 

1 

p k ' p k— 1 ] ] 

Aff k 

. fk 

Sok J 

lr_1 

L Aa k _! JJ 


At' 

2 

' V*p k+1 + 2V*P k + V*P k _! 

2 . 


L 4 


_ At (V*W k + V*W k _;L) + W 3k - ft 3k _i _ W 5k + 

- 2 2 + g (2.80) 


k= 2 (l )N 


For the top layer (k=l) its form is, 


1 1 

[ p 2 - p ll 

x °1 P NP1 * 

[At] 2 

V 2p 2 + V^Pi] 

*7 A& i 

l°2 ~ a \\ 

etT - 

2 

4 ] 


_ At (V-Wj) + W 31 _ W 51 ^ gz s 
- 2—2 MT -T- 


(2.81) 


and at the lowest layer (k=NPl) its form is. 


_1 1_ f PNPl ~ p n 1 _ P NP1 _i_ + fMl 2 T V * P NP1 + V * P N 1 

t N * A% lo NP1 - a N J RT* A<% l 2j [ 4 j 
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W 3N _ W5N At 
»N T + 2 


7-W n 


_ &Zs 


RT 


NP1 


Ld 


N 


(2.82) 


It follows from equations (2.80) to (2.82) that equation (2.75) can be 
represented by a set of NP1 Helmholtz equations. However, these 
equations are coupled and need to be made linearly independent before 
a two dimensional direct solver could be applied for their solution. 

The set of linear equations (2.80) to (2.82) can be represented 
by a single matrix equation 

AP + V 2 BP = F (2.83) 

A and B are (NPlxNPl) matrices, P is a column vector of length NP1. F 
is also a column vector of length NP1 . Equation (2.83) can be written 
as 

P + V 2 A _1 BP = A -1 F 
Let A _1 B = C 

Thus equation (2.84) becomes 

P + V 2 CP = A _1 F 

Matrix C can be diagonalized using a similarity transformation. 

If two matrices T and T -1 can be determined such that 


(2.84) 

(2.85) 

( 2 . 86 ) 


then 


where 


T C T -1 = e 


e T 


(2.87) 


( 2 . 88 ) 


e = 


C NP1 
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is a diagonal matrix representing eigen values of a matrix C. 
Substituting equation (2.88) in equation (2.86) we obtain, 

p + v 2 X _1 £ T P = A -1 F (2.89) 

Operating equation (2.89) by T, we obtain 

TP + V 2 € T P = T A -1 F 


Let T P = R and T A -1 = E (2.90) 

R + V 2 eR = EF (2.91) 

sasc met 

Equation (2.91) is a matrix representation for a set of linearly inde- 
pendent Helmholtz equations. Solution of equation (2.91) by the direct 
method provides R. Thus values of P can be obtained by using (2.90). 

Let V be the matrix of eigen-vectors of matrix C corresponding to 
eigen values e i.e. 

ass 

C V = V e (2.92) 


or 


C = V e v 


-l 


(2.93) 


But £ = X -1 £ 1 

therefore T _1 = V = Eigen vector matrix 

For finding the eigen values (e) and eigen vectors (T -1 ) of matrix C 
(=A“^B), the elements of matrices A and B are required. These can be 
written from equations (2.80), (2.81), (2.82) as follows. 

Matrix A. 

. 1 _1 1 _ 

l- 1 fT* A5 1 Aa l (2.94-a) 

A 1 1 1_ 

!.2 A a 1 Aa a (2.94-b) 

a 1,NP1 ■ Sr 

NP1 


(2.94-c) 
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V* , A»i-1 AOi 

1-1 


(2.94-d) 

111-1 

1 1 


h* A#i AOi 

1-1 

A& 2-2 Act 2 

(2.94-e) 

1 1 1 

Adi AOi 

111 


(2.94-f ) 


ANP1.NP1 = - (2.94-h) 

Matrix B 

B 1(1 = 0.25 (At/2) 2 (2.95-a) 

B 12 = 0.25 (At/2) 2 (2.95-b) 

B i,i-1 = 0.25 (At/2) 2 (2.95-c) 

Bi^ = 0.50 (At/2) 2 (2.95-d) 

B i,i+1 = 0.25 (At/2) 2 (2.95-e) 

%P1,N = 0.25 (At/2) 2 (2.95-f) 

b NP 1,NP1 =0.25 (At/2) 2 (2.95-g) 

Equation (2.91) can be written for each component, as the following 

R k (x,y,t) + e k V 2 R k (x,y,t) = G k (x,y,t) k = 1(1)NP1 (2.96) 

Here each k represents one vertical mode. It is seen that the origi- 
nal three-dimensional elliptic equation has been decoupled in the ver- 
tical and transformed into a set of two-dimensional Helmholtz 


RT^ 

NP1 


(2.94-g) 

(2.94-h) 


equations . 

The finite difference form of equation (2.96) consistent with the 
horizontal discretization is 

Rij + C V 2 R i j = Gjj i = 2(1)L-1, j = 2 (1 )M-1 (2.97) 

where 
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R. . . - 2R.. + R. - . R. . , - 2R-. + R. . , 

_ 1+1 >J i) i-l »3 i»J+l i) i,j-l 

V*r.. = + : 

ij 


A x? 


A y“ 


R i,j + i - R y-i ,an ^ 


2 Ay 


(2.98) 


For simplicity, subscript k is omitted. The values of R^ at the horizontal 
boundaries are known from the boundary conditions. 

The solution of equation (2.97) requires the following steps : 

(i) The introduction of a variable S^, defined by, 

s ij = V%- (R L ) r R i>i )(w > /<L - 1 > < 2 -"> 

This converts equation (2.97) into one with homogeneous boundary 
condition at the east-west boundaries, i.e. 


s.. + evv = H.. 
ij ij ij 

i=2(l)L-l , j=2(l)M-l 

(2.100) 

S.. = S T . =0 

1] L>J 

j=l(l)M 

(2.101) 


Here j, and Sj ^ are to be prescribed. 

and H.. = G~ - (1+eV 2 ) [Rj. - (R L j - R^-Xi-DAL-l)] (2.102) 

(ii) The variable S-. and the forcing function H.. are extended in 

xj ij 

the x-direction anti -symmetrically about the eastern boundary i=L, 
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SS 


ij 


Sij i=l(l)L-l 

' s 2L-i,j i=L(l )2L-2 


H 


HH 


ij i=l (l)L-l 


ij 


(2.103) 


(2.104) 


-H 2L-i,j 1=L(1 )2L-1 
The Helmholtz equation for SSij thus reads 

SSij + eV 2 SSij = HHij i=l ( 1 )2L-2 , j=2(l)M-l (2.105) 

SSij is periodic in the x-direction and SS-q and SSj^ are 
prescribed at the boundaries 

(iii) A discrete Fourier transform in x-direction is applied to 
equation (2.105), 

in 


1 tan4> . 

+ J- 

y 2 a Ay 


SS* . + 
J-i 


+ e 


l-4esin g 2 (L-l ) _ 2e 
Ax 2 Ay 2 


1 is' 

J J 


1 - tan(j)j 

Ay 2 2aAy 


SS j+1 = HHj 


1=1 ( 1 ) 2L-3 , j=2 (1 )M-1 
Here the Fourier transform is defined by 


FFj _ _J 

2L-2 


2L-2 r— 2I17T 
X e _y 1 2L-2 FFj j 
i=l 


(2.106) 


(2.107) 


It can be shown that the inverse transform gives 
2L-3 2J8i7T =1 

FFj j = X e/-l 2L-2 FFj , (2.108) 

(iv) Equation (2.106) represents a set of M-2 linear algebraic 

==! ==i ==i 

equations with M-2 unknowns, i.e., SS 2 , SS 3 , .... SSjj-i> which can be 

solved by standard inversion methods. In fact, the matrix of equation 

(2.106) is tridiagonal, and is solved by a Forward-Backward iteration. 

==l 

(v) After solving SSj from equation (2.106), one obtains SS^ from 
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equation (2.108). 

(vi) Next we calculate Rjj from equation (2.99). 

(vii) Finally we calculate Pjj^ from equation (2.90). 

We shall next address the Lagrangian interpolation. 

The semi -Lagrangian advective scheme was first proposed by 
Krishnamurti (1962). A linear stability analysis of this scheme was 
carried out by Mathur (1970). A multi-level formulation of the semi- 
Lagrangian advection scheme explored the formation of an intertropical 
convergence zone, Krishnamurti (1969). In these studies an Euler- 
backward time differencing scheme was used. An excellent review of 
the aforementioned efforts was presented by Haltiner (1971). Rapid 
further developments on the formulation of semi-Lagrangian - semi- 
implicit versions of primitive equation models emerged from a series 
of recent studies, Robert et al. (1985), Ritchie (1986), Bates and 
McDonald (1982). In these studies the use of semi-implicit time dif- 
ferencing and long local trajectories for the horizontal advection 
resulted in a drastic reduction of the computation time. Proposals 
for the removal of the Helmholtz equation for the pressure via a 
transformation of variables was made by Tanguay and Robert (1986). 
This can, in principle, provide a further reduction in the com- 
putational time. A local analytical representation of the basic 
variables and forcing functions is normally done via use of local 
interpolation surfaces. Ritchie (1986) proposed yet another improve- 
ment in the mathematical formulation of the semi-Lagrangian advection, 
that permits the removal of such an interpolation phase. A proposal 
to cast this as a spectral-semi-Lagrangian - semi-implicit problem was 
recently suggested by Ritchie (1987). Many of these newer methods 
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remain to be tested at the present time for possible implementation in 
multi-level models. Mathur (1983) proposed a one-step semi-Lagrangian 
model that utilizes an explicit time differencing scheme. The present 
study utilizes a one-step semi-Lagrangian advection with a semi- 
implicit formulation for the time differencing scheme. 

For any variable A, the general time dependent problem can be 
written as: 


dA _ p (2.109) 

dt 


where F is the forcing function and 

4- , i. «. u 8_ + v a_ 

dt 3t 3x 3y 

The numerical integration of equation (2.109) by the semi- 
Lagrangian advective scheme is carried out in the following manner: 

(a) over the time At the origin of a parcel reaching a certain 
grid point is traced. 

(b) The value of A at the new location is advected from its place 
of origin. 

(c) an average forcing F is obtained over the period At and is 
added to the advected value of A. 

Mathematically this is expressed by 

Ag(t+At) = Ap(t) + FAt (2.110) 

i.e. the air parcel starting at point P in time At reaches the point Q 
and the average forcing in the time At is F. 


where 


t+At 



Fdt 


( 2 . 111 ) 


To numerically integrate equation (2.109) in order to obtain 
equation (2.110) certain approximations have to be made. Given the 
velocity fields at the time t and assuming that the forcing over the 
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period t to t+At is constant the distance travelled in time At is 
given by 

Ax = u(t)At + 3gF u ( At ) 2 (2.112) 

Ay = v(t)At + %F v (At) z 

Where Ax and Ay are the distances travelled along the x and y 
directions respectively. This approximately gives us the location of 
the point P. Using an interpolation technique the value Ap(t) is ob- 
tained. To obtain an average value of F (equation (2.111)) we assume 
that this forcing is also constant and is equal to the forcing at the 
parcel at point P. i.e. 

F = Fp(t) (2.113) 

A second order Lagrangian interpolation scheme is used in this 
model. This scheme uses a nine point stencil. Higher order schemes, 
at the expense of considerable increase in computational time could be 
designed. The accuracy of a given scheme depends upon the order of 
interpolation used, both to find position of, and the value of the 
fields, at the departure points. McDonald (1987, 1984) has analyzed 
the accuracy associated with the different spatial interpolation sche- 
mes. More recently McDonald and Bates (1987), Temperton and 
Staniforth (1987) have suggested the way to improve the estimate of 
the departure point positions. In general these schemes require 
information from more than two time levels. 

In the x-momentum equation the grid point values of the meri- 
dional velocity v are staggered with respect to the points that carry 
the zonal velocity u, hence certain averaging procedures are required. 
The velocity component v at grid point is obtained as follows: 

^ijk = %< v i-ljk + v i jk + v i-lj+ik + v ij+lk> (2.114) 
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and a forcing function corresponding to the v component is averaged 
from: 


F ijk = ^(Fi-ljk + F ijk + F i-lj + lk + Fij + lk) (2.115) 

u 

Since the U£j k and Fjj k are known the distance travelled in time 
At are given by: 

u 

Ax = Uij^At + ^F 1 j k (At) 2 (2.116) 


Ay = v ljk At + 3gFi Jk (At)* (2.117) 


This gives us the location of point P for the semi-Lagrangian 
form of x-momentum . 

For y-momentum equation since the grid points with values of u 
are not located at the V points, certain averaging is again required. 
The velocity component u at Vjj k grid point is obtained as follows: 

Uijk = %( u ij~lk + u i+lj-lk + u ijk + u i+ljk) (2.118) 


and the forcing function corresponding to the u component is averaged 
by: 

Fijk = ^(Fij-ik + F i+i j-ik + F ijk + F i+ljk) (2.119) 

again the distance travelled are given by: 

Ax = ujj k At + % F l j k (At) * (2.120) 

Ay = v ijk At + % F j jk (At ) 2 (2.121) 

This gives us the location of point P for the semi-Lagrangian 
form of y-momentum equation. 

We shall next address the thermodynamic, continuity and moisture 


conservation equations. 



Since the temperature, surface pressure and the moisture 
variables are carried at the same grid point a general procedure for 
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all three equation is given here. 

Again since for these variables the u and v do not fall on their 
grid points, the following averaging procedures are employed. 


Q ijk = ^( u ijk + u i+ljk) (2.122) 

F ijk = ^( F ijk + F i+ljk) (2.123) 

^ijk = %( v ijk + v ij+lk) (2.124) 

F ijk = ^( F ijk + F ij+lk) (2.125) 

Knowing these values of above variables the distance travelled 
can be obtained as: 

~u 

Ax = u ijk At + %F ijk (At) 2 (2.126) 

Ay = v iJk At + %Fij k (At) B (2.127) 


This gives the location of point P, for the semi-Lagrangian form 
for the respective equations. 

The integration involves the following steps: 

1. The forcing, sj are calculated at grid points and are stored in S j . 

2. Ax and Ay are calculated using the procedure discussed above. 

3. Forcings s j , at the grid points, are interpolated to the origin of 
the parcel using Ax, Ay by a similar interpolation scheme and stored 
in Sj °. 

4. Forcings for the elliptic equation are calculated using Xj's and 
stored in Q. 

5. The three dimensional elliptic equation is transformed into a set 


of 2 dimensional linearly independent Helmholtz equations and these 
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are solved by direct method to obtain P + . 

6. Given P + and W + , the prognostic and diagnostic variables are calcu- 
lated as follows: 


u-component of wind 


w + _ Ai 3£ + 

1 2 3x 


v-component of wind 
v+ = W 


r+ _ w+ _ At 3P + 


2 2 3y 

Surface pressure 

+ l + 

lnp s = ^T-*(P S - g z s) 

S 

Moisture Variable 

q + = W + 

4 


Geopotential Height 

Z+ = g (P+ ~ RT * lnp s) 

Temperature 

T + _ _SE az+ 

R 3a 

Vertical velocity a 

° = it n w s - inp s + - ^ v *^ + ] ** 

(7t 

These steps are repeated to obtain the forecast. 


(2.128) 


(2.129) 

(2.130) 


(2.131) 


(2.132) 


(2.133) 


(2.134) 


2.4 Initialization : 

Sugi (1986) proposed a dynamic normal mode initialization 
procedure that utilizes a forward backward iteration (around the first 
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time step). In this method the forward-backward iterations are 
carried out with the linear terms, while keeping the physical and 
nonlinear effects fixed. After every 100 iterations the physical and 
nonlinear terms are updated and this process is continued until 
roughly 5 iterations are completed. Sugi has demonstrated a rapid 
damping of high frequency modes with a solution which corresponds clo- 
sely to that obtained from nonlinear normal mode initialization with 
physics. Sugi also demonstrated a rapid convergence in the time 
history of mean square surface pressure, mean square divergence and 
grid point values of surface pressure. Sugi's study was carried out 
within a global spectral model at a resolution of 42 waves 
(triangular). The transform grid at the resolution corresponds 
approximately to 250 km 2 . When the Sugi iteration is carried out in 
regional grid point models at a resolution of roughly 50 km 2 or higher 
a rapid damping of the gravitational high frequency modes occurs. 
However in mesoscale models with still higher resolutions this itera- 
tion scheme appears to fail and requires a reformulation. The present 
study utilizes grid sizes of the order of 50 km or higher and this 
method appears to converge quite rapidly. 


2.4.1 Formulation : 


The linear and nonlinear terms in the model 
separated for the Sugi initialization as shown below: 
(a) X-momentum equation 


Du 

Dt 


3u 

a ~ + v 
8a 


f + - tancj> 
a 


■j- 


RT' |^1 np s + Fx 


equations are 


3P (2.135) 
dx 
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where D = 3 + u3 + y 3 

Dt 3t 3x 3y 

equation (2.135) can be rewritten as: 

3u 3P 


(2.136) 


+ vf + 


3t 


3x 


a - u |S - v JJ + 2S tan<J> - RT' fjlnp, + Fx 


3a 


Bx 


= F u + F u 


i n 

(2.137) 

I n 

where F u and F u are the linear and nonlinear parts of the equation 


respectively and are given by: 


u 

n 

? u 


3P 

3x 


+ vf 


(2.138) 


-a^-u|^-v-|^ + — tan(J» - RT' lnp s + Fx , . 

3a 3x 3y a 3x (2.139) 


(b) Y-momentum equation 


Dv 

Dt 


‘ 3v 
3a u 


f + - tan<J> 
a 


- RT ' fe ln Ps + F Y - 


3P (2.140) 

3y 


following the same procedure we can write equation (2.140) as: 

- v If - ^ tani » - RT ’ f^ ln p s + p y 


3v 

BP 

- uf + 


u |2 
Bx 

at 

ay 

- (T — - 

Bcr 

or 


3v _ F 1 

n 




3t Fv 

+ F V 


where 






i 

?v 

n 


(2.141) 

(2.142) 

(2.143) 


3P 
3x 

Fv «- a |2. u |v_ v |v_uu _ |_j „ F „ 

v 3a 3x 3y a 3y s Y (2.144) 

are the linear and nonlinear terms respectively. 

(c) Continuity Equation : 

In the Eulerian form, the continuity equations can be written as: 


3t 


+ v»V lnp s = 


3d 

3a 


- V»v 


(2.145) 


where 


v = iu + jv 


(2.146) 
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and 


V = i f + j I— 
3x J 3y 


(2.147) 

Integrating equation (2.145) with respect to a from <7=1 to a=a T we get: 

1 1 1.1 
| -|^lnp s da + | v*Vlnp s da = - j ~ da - j V*v dcr 


CT'p CT^p 

4 


CTj 


0T 


defining v = 


( 1-C7j) 


a T ) } 


\fdcr 


ax 


Using the boundary conditions j = CT 1 - = 0 

(Jj 


we can write equation (2.148) in the form 
B 

— lnp s = - v»Vlnp s - V*v 

or 


3 

at 


lnp s = F 1 


+ f ds 
ps P s 


where 


(2.148) 


(2.149) 


(2.150) 


(2.151) 


ps 


V*v 


F ps = “ S*Vlnp s 

are the linear and nonlinear terms respectively. 

Subtracting equations (2.145) and (2.150) we can derive a 
diagnostic equation giving 


(2.152) 

(2.153) 


aa 

acr 


= (v - v)»Vlnp s + V«(v - v) 


(2.154) 


(d) Thermodynamic equation : 


Dt [ 7" H + <Tlnp = ] + c ^ 


KRT' 
*r 

cry 


3a 

3a 


+ V«v 


+ Q 

(2.155) 
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1 0P 

let X = —y — we can write (2.155) as: 
y 8a 


DX D 


+ a — 3np,, + a = - 


Dt v Dt 


kRT’ 


y cry 


aa 

3ct 


+ V* v 


+ Q 


(2.156) 


substituting for |-£lnp s from equation (2.150): 


DX 

Dt 


+ a 


|^lnp s + v.Vlnp s ] + a = - ^ [ §§ + V»v ] + Q 


(2.157) 


or 


DX 

Dt 


+ g(- $*Vlnp s - V»v + v*Vlnp s ) + a 


y 'g + icRT ' 
7 cry 3 * 


8a 

8a 


+ V»v + Q 


(2.158) 


or DX 
Dt 


= aV*v + 


p(v - v) *Vlnp„ - ^ KRT ' 


+ sr 

y gy 


3cr 

8g 


+ V*v + Q 


- g 
2.159) 


or DX 
Dt 


= gV*v - g + 


r(v - v)»Vlnp s - 


g kRT 

■» + s 

y cry 


3g 

3g 


+ V»v 


+ Q 

(2.160) 


or = F* + F n 
Dt X X 


(2.161) 


where 


F x = gV*v - g 


(2.162) 


F n = - g ( v - v)*Vlnp s _ Xlg + 
x y gy 


3g 

3g 


+ V»v 


(2.163) 


are the linear and nonlinear terms respectively. 

(e) Moisture conservation equation : 

Since the dynamic initialization procedure does not effect the 
moisture variable in a direct way, no attempt is made to separate the 
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moisture equation into linear and nonlinear parts. 

2.4.2 Procedure: 


First the model is run in the complete form and updated fields 
are obtained at first time step. From the fields stored at the time 
zero, the total tendencies are obtained. The nonlinear tendencies can 
be readily obtained by subtracting the linear tendency from the total 
tendency. Following Sugi (1986) these nonlinear tendencies are held 
fixed through the forward-backward integration. This forward-backward 
integration is performed for n times, after which the updated initial 
fields are used to run the model for one time step. This rerun of the 
model gives a new set of nonlinear tendencies, which are again held 
fixed through n forward-backward integrations. This procedure is 
repeated for several times. Fig. 3 gives a flow diagram for the 
entire procedure. 

The linear integration scheme is different from that used in the 
model. In the model a semi-Lagrangian, semi-implicit scheme is used 
for the spatial and temporal time differencing. The linear forward- 
backward steps are performed in the Eulerian sense using the selective 
damping scheme proposed by Okamura (see Haltiner and Williams page 
387). The Okamura scheme consists of a forward followed by a backward 
step and averaging. 

F n+1 = F n + At [ at ] n _ (2.164) 


new 1 

Fn = Fn + 1 - At 



(2.165) 
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Figure 3: Flow diagram of the initialization procedure. 
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F n = 3F n - 2Ffj V (2.166) 

where F maybe u, v, lnp s or x and n stands for the time step. 

The variable x is calculated in the very beginning of the proce- 
dure from its definition. 

2.5 Boundary Conditions : 

a) Vertical : 

In the vertical o is set to zero at model top and bottom for all 
the times t. The lowest level of model follows the surface terrain. 
These conditions can be represented as: 


at 

r-t 

O 

II 

b 

and 

a = 1.0 

and at 

a = 1.0 


Z (x ,y ,a=l ) = Z s (x,y) 

where Z s is the surface terrain field and is function of x and y. 
b) Horizontal : 

Numerical treatment of lateral boundary is a difficult but a very 
important aspect of limited area modelling. The flow in the interior 
of the domain is sensitive to the specification of the lateral boun- 
dary conditions, even for the short range forecasts. 

An excellent review of commonly used lateral boundary conditions 
for the limited area models is given by Davis (1983), Sundstorm and 
Elvins (1979). The necessarily limited domain introduces the lateral 
boundary conditions within the atmosphere whose physical interpretation 
is obscure. Further incorrect specification of lateral boundary 
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conditions itself can generate local unstable modes even if the basic 
finite difference approximations used in the interior are stable. The 
uniqueness and well-posedness of lateral boundary conditions of a 
limited area model are still open questions. However over the years 
many pragmatic approaches have been successfully developed and applied 
towards the short range forecasting. These could be broadly put in 
the two categories. The first approach is so called grid nesting 
where a much larger domain is also modelled using a coarser resolu- 
tion, and the evolving fields are then used to determine the boundary 
conditions for the smaller domain (Chen and Miyakoda 1974, Miyakoda 
and Rosati 1977, Leslie et al. 1981). The second approach is to for- 
mulate open boundary conditions which allow the disturbances to 
radiate out of the modelled domain with minimum boundary reflections. 
The most used procedures under this category are diffusive damping 
schemes in which the boundary noise problem is dealt by introducing a 
marginal zone of large diffusion (Mesinger, 1977), tendency modifica- 
tion schemes in which the time tendencies are assigned a weighted 
average of the externally specified fields and internally forecasted 
fields (Perkey and Kreitzberg 1976), relaxation or nudging schemes in 
which the flow near the boundary is relaxed towards a large scale ana- 
lysis through the use of Newtonian or diffusion term (Kallberg and 
Gibson 1977, Tatsumi 1980). Apart from these many explicit radiation 
schemes based on outflow-inflow conditions have also been used 
(Williamson and Browning 1974, Okamura 1975, Orlanski 1976, Miller and 
Thorpe 1981, Kurihara 1983). 

A similar pragmatic approach has been used in the present model. 
Two such options of the boundary conditions are currently available: 
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I. Fixed boundary conditions : In this case all the prognostic 
variables i.e u, v, p, lnps and q are held fixed in time. To avoid 
the wave reflection, caused by such a boundary condition, a Laplacian 
type smoother near the boundary is applied. The same smoother is also 
applied in the interior of the domain to suppress the noise. However, 
as the lateral boundaries are approached, the smoothing coefficient 
increases. For the 1.875 lat/long grid the smoothing coefficient 
nearest to the boundary is 16 times its value in the interior of the 
domain. This increase is a function of grid resolution. For smaller 
grid resolutions this multiplying factor has a smaller value. 

II. Time dependent boundary conditions : Here the prognostic 
variables are allowed to change with time. The future values of the 
prognostic variables can be obtained either from a global forecast or 
from the observations. These specified values then are merged with 
the forecast produced by the model. In the present study we use time 
dependent boundary conditions based on future observations. If at any 
time t the model forecast is Xf and the externally specified value is 
x s , where x is any prognostic variable, then the following expression 
is used to give the merged forecast 

x c = (1-a) Xf + ox s 

where a is a constant and equal to 1.0 at the boundary. Away 
from lateral boundary it decreases and becomes zero in the interior. 
The specific values for a are 1.0, 0.45339, 0.20556, 0.09320, 0.01916, 
and 0.0 for grid points, 0, 1, 2, 3, 4, 5 and more, distance away 
from the boundary. 

The above procedure is carried out for u, v, z, lnps and q on the 
a surfaces and is repeated at each time step of the model 
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integration. The specified value X s is obtained through a linear 
interpolation between two observation times. 



CHAPTER 3 


PARAMETERIZATION OF PHYSICAL PROCESSES 


Apart from the dynamical aspects, the model also includes a 
number of physical processes. These comprise of: 

• Large scale condensation, 

• Shallow moist convection, 

• Deep cumulus convection, 

• Planetary boundary layer, 

a) Surface fluxes 

b) Fluxes above the surface layer. 

• Radiative Processes 

a) Short wave radiative processes, 

b) Long wave radiative processes, 

• Surface energy balance, 

• Orography . 

A fairly accurate parameterization of the above mentioned physi- 
cal processes in an NWP model is very essential to capture the various 
details of importance in a forecast. Following is a brief description 
of the parameterization schemes used in the model for all individual 
processes. Numerical algorithms are written in a modular form for 
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each of these physical processes. These parameterizations provide the 
physical forcing terms for the previously discussed equations of 
motion. 

3.1 Parameterization of Large Scale Condensation : 

Large scale condensation within the model is essentially the 
removal of supersaturation at any grid point of the model. This 
excess moisture develops due to the other dynamical or physical pro- 
cesses such as vertical advection of moisture fluxes and radiative 
cooling, etc. 

Super saturation at each time step at a given grid point is 
defined by 

Aq = q-q s ( s.i ) 

where q and q s are the specific humidity and saturation specific humi- 
dity at that point. 

If Aq>0, then at that le 1 of atmosphere the excess humidity is 
precipitated out. The heating due to the precipitation ( LA< 3/ At ) is 
computed and added as a source term in the thermodynamic equation. 

c p 1 M = • . . • + Ma 

e at At (3.2) 

Similarly change in moisture due to supersaturation is also added to 
the moisture continuity equation. 

dq =••••- Aq 

3t At (3.3) 

Thus the supersaturation is simply condensed out with an equiva- 
lent heat released at that level of the atmosphere. 
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Calculation of the saturation specific humidity q s is usually 
carried out from the use of various approximations such as the Tetens 
formulae, q s (T) = 0.622 e s (T) (3.4) 


0.622 e s (T) 
p - 0.378 e s (T) 


vapor pressure 


6.11 exp 


a (T-273, 16) 
(T-b) 


(3.5) 


where the constants a and b in the expression for vapor pressure are 
defined in terms of saturation over water (a = 17.26, b = 35.86) or 
over ice (a = 21.87, b = 7.66). The Tetens formulae have been tested 
and found to be reasonable for the construction of moist adiabats in 
the troposphere. 


3 . 2 Parameterization of Shallow Convection: 


Following Tiedke (1984) we assume a simple diffusive form for the 
parameterization of shallow convection. Shallow convection is invoked 
over conditionally unstable undisturbed situation (rainfall rate < 5 
mm/day) via a K theory. The humidity and the thermal equations are 
expressed by. 


3q = .....3 K n 9q 

at ap Q dp (3.6) 

and 

38 = 3_ K e 96 (3.7) 

at ap 3p 


Furthermore, following Tiedke (1984), we set, Kq = = Kg 2 p 2 , 
the coefficient varies with height as the square of the density. 
Based on numerous studies with the GATE, ATEX and AMTEX data sets 
Tiedke suggests a value of K = 25 m 2 s -1 . However, setting a constant 
value of K does not provide a smooth continuity with the fluxes at the 
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lifting condensation level. It is desirable that Kg^p^ at the 

3p 

lifting condensation level be consistent with the planetary boundary 
layer flux at that level. If we define: 


K = F q (LCL)/L 
p 2 g is 

0 p)lcl 


(3.8) 


then the humidity flux in the subcloud layer and the shallow convec- 
tive cloud layer are smoothly connected. A similar definition of K 0 
is also used for the same reasons. Above the lifting condensation 

level in the cloud layer K 0 it is negative, implying a downward eddy 

9p 

flux of sensible heat from the base of inversion to the cloud base. 

On the other hand 3cj; essentially positive in the shallow convec- 

3p 

tive cloud layer. Since the humidity flux is small both at the cloud 

base and at the cloud top, over the lower half of the clouds — K 8q 

3p 14 3p 

is negative and is positive over the upper half of the clouds. Thus 
these shallow convective clouds tend to moisten the layers of the 
atmosphere over the upper half of the clouds. Heckley (1984) has sum- 
marized the results of recent numerical predictions made with the 
ECMWF model . The experience at ECMWF shows that the inclusion of 
shallow convection improves the overall tropical forecasts slightly. 
In its absence the planetary boundary layer tends to be too moist and 
the region below the inversion layer tends to be too dry. Although 
the present formulation is very simple it is shown to provide an 
overall improvement of the thermodynamical structure over the 
undisturbed areas. In the present formulation the vertical eddy flux 
of, momentum is included in the planetary boundary layer; however, it 
is not included within the cloud layer. 
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It is easy to see that moisture is conserved by this parameteri- 
zation of shallow convection. The total convergence of flux of humi- 
dity in the cloud layer is equal to Fq(LCL), i.e. the flux at the 
cloud base. 

It is important to assure that these shallow clouds are nonpreci- 
pitating and do not systematically produce saturation on the large 
scale. To assure that it is desirable to set an upper limit on the 
final specific humidity q F produced by the process, 
i.e. qp/q s * 08 

This is easy to implement formally although it violates the conser- 
vation slightly. However, this is necessary in order to avoid the 
eventual formation of layer clouds and associated large scale conden- 
sation. 

The present formulation of shallow convection is very simplistic, 
and deserves to be investigated in further detail. 

3.3 Parameterization of deep moist convection : 

The procedure followed here is identical to that used in 
Krishnamurti et aK (1983a). Deep convective parameterization is 
invoked if the atmosphere is conditionally unstable and if a net 
supply of moisture convergence is available. The existence of con- 
ditional instability is determined from vertical gradient of moist 
static energy i.e., 

§- (gz + CpT + Lq) < 0 


(3.9) 




50 


The total supply I may likewise be split into the moistening and 
heating parts by the relation, 

Iq = lb = I L b(l+rj) (3.16) 

and Ifl = I(l-b) = I L (l-b) (1+n) (3.17) 


The thermal and the humidity equations are expressed by: 

M + VV0 + = a fl ( e s~ e + (3.18) 

9t 3p a t 3p 

IS + v*Vq = a q f q S -q ) (3.19) 

8t AT 

Where ag and a q are defined by the relations, 

a e = £e = I(l-b) = I L (l+n)(l-b) (3.20) 

Qe q 9 Qe 


»q “ Is = lb = lLb(l +i?) 

Qq Qq Qq 


(3.21) 


The parameterization is closed if b and 17 are somehow determined then 

ag and a q may be evaluated from the relations, 

a e « l£ - LU-kl = I L (l-b)(l+n) (3.22) 

Qe Q 0 Qe 

(3.23) 


and 


a q = = lb = i L bd+n) 


Qq Q Q 


Krishnamurti et ad. (1983a) proposed a closure for b and rj based 

on a screening multiregression analysis of GATE observations. 

— = ajf + bjw + cj (3.24) 

J L 

E_ = a2f + b2<«) + C2 (3.25) 

J L 

Where aj , bj , cj , a2> b2> and C2 are regression constants whose magni- 
tudes may be found in Krishnamurti et al. (1983a). C and u are respec- 
tively the relative vorticity and the vertically integrated vertical 


velocity. 



51 


Thus, In numerical weather prediction, £ and « determine M/Il and 


M_ = b ( 1 +rj ) (3.26) 

IL 

B_ = ( l+T}) ( 1-b) (3.27) 

It 

These two relations determine b and r ?. 

We can find ag and a q from the relations, 

_ l L (l-b)(l+r}) _ R_ 
a e ~ ~ 

Qe Q * 

and 

. _ It b(l+7)) _ M_ 

a q - - 

Qq Q( 3 

It is also of interest to note that the apparent heat source Qj and 
the apparent moisture sink Q 2 , for this formulation, may be expressed 
by: 

Ql = a e f c P 1 ( e s-6) + (tfCp 1 Ml + C p 1 [Hr + H s ] < 330 > 

[ 6 At 0 3pJ 0 

q 2 = -l F a q (qs~Q) + w l3l (3.31) 

l At 3pJ 

The code monitors the apparent heating and moistening for output 

CENH'L u'RARY 

1 * 


(3.28) 


(3.29) 


R/Il; 

since 


and 
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purposes. The total convection precipitation is given by, 

PB 

P T = I / Cp T (e s -e) dp (3.32) 

g Px 6 At 

3.4 Parameterization of surface fluxes 


Surface fluxes in the model are calculated following the for- 
mulation of Chang (1978). The formulation is based upon the assump- 
tion of the existence of a constant flux layer of about 50 meters 
depth near the earth surface in which the fluxes of momentum, sensible 
heat and latent heat (moisture flux) are constant with height. 

The surface fluxes of momentum, heat and moisture are given as 



F m = -p<u'w’> s 

(3.33) 


Fs = -pc p <e'w'> s 

(3.34) 

and 

F l = -pL<q'w' > s 

(3.35) 


where u', w' , 0' and q' are the deviations of horizontal velocity, 

vertical velocity, potential temperature and specific humidity from 
their large scale values. Thus these are associated with the tur- 
bulent eddies. The < > is a large scale time average. Subscript s 
denotes the surface layer. The computation of these surface fluxes is 
based on similarity theory described by Chang (1978). 

The basis of the similarity analysis are the planetary boundary 
layer observations, e.g. Businger, et al. (1971). According to these 
observations, nondimensionalized vertical gradients of large-scale 
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quantities such as the wind (u) , potential temperature (0) and 
specific humidity (q) in the surface layer can be expressed as univer- 
sal functions of a nondimensional height (z/L) where z is the height 
above the earth's surface and L is the Monin-Obukhov length. 

Monin-Obukov length is defined by, 

L = u*Vk/S0 # (3.36) 

Here 0 = £/0 o where 0 O is the mean potential temperature in the 
constant flux layer. u* and 0* are related to the surface momentum 
flux and surface sensible heat flux by 

u z # = - <u'w'> s (3.37) 

and u*0* = - <0'w'> s (3.38) 

We express the nondimensional wind shear, potential temperature gra- 
dient and moisture gradient by the relations, 



kz 3u 

= <f) m (z/L) 

(3.39) 


u* 3z 




kz 80 
0* 3z 

= <|> h (z/L) 

(3.40) 

and 

kz 3q 
q* 3z 

= <J> q (z/L) 

(3.41) 


The empirical fits of the boundary layer observations are usually 
separated in terms of stability. Stability is expressed in terms of 
the sign of the Monin-Obukhov length L or that of the Bulk Richardson 
Number Rig. 

d0 

R _ « dz (3-42) 

RiB ~ P 

du] * 
dzj 

For the stable surface layer 


L > 0 


(3.43) 
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or R iB > 0 (3.44) 
For the unstable surface layer 

L < 0 (3.45) 
or R 1b < 0 (3.46) 
Since, apriori the Monin Obukhov length is an unknown quantity, sta- 
bility is assessed from the sign of the Bulk Richardson Number Rj B . 


Empirical fits on the boundary layer observations (Businger et 
al . , 1971) for wind shear and potential temperature gradient for 

stable and unstable cases, are expressed, following Chang (1978), by 
the following nondimensional relations, 
a) Stable case 

kz/u* |~ = 1-0 + 4.7 z/L (3.47) 


kz/e* ~ = 0.74 + 4.7 z/L 

and 

b) Unstable case 


(3.48) 


kz/u* |S o (i - « z/L) -3 * (3.49) 

0 z 

kz/e* §£ = 0.74(1 - 9 z/L)- 3 * (3.50) 

u# and 8* within the constant flux layer are constant with height 
since fluxes are constant. Using this fact and the definition of 
Monin-Obukov length given by (3.36) we can integrate the equations 
(3.47) and (3.48) to provide u* and 6* for stable case as: 


k(u 2 - uj) 
lntzg/zj) + 4.7~Az 
L 


(3.51) 
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and 


t> _ k(e 2 - Oj) 

* ~ 0.74 inlz^/zj) + 4.7 Az (3.52) 

L 

where , u 2 and 0^ , are wind speeds and potential temperature at 
heights Zj, z 2 in the constant flux layer. Az is equal to (z 2 - Zj). 
Using (3.36), (3.51) and (3.52) by eliminating u* and 6* a quadratic 
equation in Az/L is obtained. Its one solution which satisfies the 
physical constraint that L>0 for the stable condition is 


Az = In 
L 


(Z 2 ) [ 9.4 R iB - 0.74 + (4.88 R iB + 0.5476)% 
zj [ ST* - '44 .18 R iB . 


(3.53) 


where Rj B (Bulk Richardson Number) 

= 0 (0 2 - 0j)Az 

(u 2 - Uj) a 

Equation (3.53) also requires that 0 < Rj B < 0.212 to ensure Az/L > 0 
for R iB > 0. For stable case Az/L is calculated from (3.53) and this 
in turn provides u* and 0* using (3.51) and (3.52). u* and 0* pro- 
vide the momentum flux and sensible heat flux through the relations 


(3.37) and (3.38). 

For unstable case the integration of (3.49) and (3.50) provides: 


u, = k(u 2 ~ Uj) (3.54) 

1 R ( z 2 ) ~^1 (^2 ~ z 2 ) 

Zj ~T Zj 

and 

6* = . M0 2 - ®i) 

0.74 T in(Z 2 ) - W 2 ( z 2 > z 2 ) (3.55) 

zj TT zj 


The functions iJ/j and <|f 2 are given by 
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*M Z 2 • £2) 

IT - Zj 


In 


- 2 


| (1 + E 2 ) J 2 j (1 + E 2 Z ) 

1 rr t e 1} 1 1 n + ejt * 

tan _1 E 2 + 2 tan^Ej 


»] 


(3.56) 


and 


where 


<J>2 (£2 > £2.) 

L Zj 


2 In (1 + T 2 ) 
1 + *1 


El 


e 2 


*1 

*2 


f 1 

(! 

(1 


- 15Z 2/ L) 

- 15 Zg)^ 

9z 2 /L 

z^7zj) 

(1 " 9 £2) 

L 




% 

% 


(3.57) 

(3.58) 

(3.59) 

(3.60) 

(3.61) 


Using (3.36), (3.54) and (3.55) we can get an expression for z 2 /L. 
Chang (1978) has given a polynomial least square fit to obtain z 2 /L 
from the modified bulk Richardson number Ri m defined by 


R im = R iB z 2 
Sz 


(3.62) 


The equation is given by 

In f-z 2 ) = a 0 + a 1 ln(-R im ) + a 2 [ln(-R im )] z + a 3 [ln(-R im )] 3 
T 

(3.63) 


The regression coefficients a 0 , a 1 , a 2 and a 3 are functions of 
ln(z 2 /z 1 ) . 

Thus, knowing z 2 /L, u*, 6*, surface flux of momentum and surface 
flux of sensible heat can be obtained. 



The fluxes can also be expressed in terms of bulk aerodynamic 
formulae i.e. 


f M = PCp ( u 2 ~u 1 ) 2 (3.64) 

F S = pCpC H (u 2 -u 1 )(e 2 -e 1 ) (3.65) 

F l = pLC Q (u 2 -u 1 )(q 2 -q 1 )*GW (3.66) 

where Cp , Cy and Cq are exchange coefficients for momentum, heat and 
moisture and GW is a ground wetness parameter. By comparing equations 
(3.64), (3.65) with equations (3.37), (3.38) using the definition of 

fluxes from (3.33) and (3.34) the expressions for the stability depen- 
dent exchange coefficients can be written down as: 


a) Stable Case: L > 0, 

Momentum: 


C D = 


k® 


In 

’ z 2 ] 

+ 4.7 Az 


L z iJ 

| L 


(3.67) 


Heat: 


C H = 


-k 2 


In ^2 + 4.7 M 


0 . 741n £1 + 0 .47 M 

Z 1 L 


Z 1 L 


(3.68) 


Moisture : 

The exchange coefficient for moisture C q is simply defined by: 

Cq = 1.7 C H (3.69) 

b) unstable case: L < 0, 

The exchange coefficients are 


Momentum : 
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Heat: 


k 2 


Ch 0.74 [in z 2 /zj - iJx-jJ [in Z 2 /Z 1 - 


(3.71) 


Moisture : 


Cq = 1.7 C H (3.72) 

Ground wetness parameter GW required for the calculation of moisture 
surface flux is defined by the following emperical relations : 


GW 


a, 


■max 


- a 


a, 


max 


*min 


(1 - RH) 


(3.73) 


or GW = 0.85 (1.0 - exp(-200 . 0(0 . 25-a) 2 ) ) (3.74) 

Here RH denotes the relative humidity of air at level 2. a^x, a 
and otjnin respectively denote the maximum surface albedo of the domain, 
local surface albedo and the minimum surface albedo of the domain. 
Over ocean GW is set to 1.0. 

In practice the height zj is taken to be the roughness height zq. 

Over the oceans we use Charnock ' s formula to define roughness 
height, i.e. 

z 0 = 0.04 * u* 2 /g (3.75) 

Since z 0 is a function of the momentum flux its value presupposes a 
knowledge of u*. u* is first calculated using constant exchange coef- 
ficient in the Bulk Aerodynamic formula. This provides a first guess 
value for z Q . Later it is updated by iterating (3.75) with the stabi- 
lity dependent calculation for u* in the similarity theory. 

Over land areas the roughness height z 0 is determined using the 
formula given by Delsol et al . , 1971. 

z 0 = 0.15 + 0.2 (236.8 + 18.42*h) 2 * 10“ 8 m. (3.76) 
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Here h denotes the terrain height in meters. In our calculations 
we have set an upper limit for z 0 =l meter. 


3 • 5 Parameterization of fluxes above the surface layer : 


Above the surface layer the fluxes of heat, moisture and momentum 
are calculated following the K-theory (Louis, 1979) where the eddy 
diffusion coefficient K depends on the mixing length, the vertical 
wind shear and the stability of the atmosphere. 

The convergence of flux at any level is calculated from a diffusion 
equation, 

12 = I 3_ 

Bt p Bz 

where Q stands for u, v, B or q. At the upper boundary flux is 
assumed to be zero and at the lower boundary it is equal to the sur- 
face fluxes as calculated in section 3.4. We define the eddy dif- 
fusion coefficients for heat (Kg ) , moisture (Kq) and momentum (K m ) by 
the relations, 


pK q fi 


(3.77) 



Ke = Kq = 1* 

BV 

8z 

Si (Rj ) 

[r] 

and 

K m = ** \ 

3Vl 
Bz | 

s 2 ( R i) 

[r] 


(3.78) 

(3.79) 


! = mixing length 

Sj and S 2 express the dependence on Richardson number Rj , and are 
expressed by the following emperical relations : 

1.0 


®i- s 2 = (I rs Rp ? for Ri *° 

0 _ ( 1.0 + 1 . 286 |Hi 1^ - 8Rj ) 

b l 

(1 + 1 . 286 ^ 1 %) 


(heat, moisture, and (3.80) 
momentum ) 


Rj<0 (heat, moisture) (3.81 ) 
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and S = -° + 1 1 T4£]Rjjj ~ 8Rj) 


(1 + 1.746|R i j^) 


R i = 


ae v 

^sr 


Rj<0 (momentum) (3.82) 


(3.83) 


av 


3z 


where 0 V is virtual potential temperature. 

! denotes the mixing length and is defined (following Blackadar, 1962) 

by 


kz ; i inf = 150m (momentum) 

+ kz = 450m (heat, moisture) (3.84) 

*inf k = 0.35 Von Karman constant 


3 .6 Shortwave Radiation 

Solar radiation absorbed at the earth's surface and in the 
atmosphere is the initial source of energy causing atmospheric 
motions. A reliable treatment of solar radiation in numerical models 
of the atmospheric circulation is required. 

The main absorbers of the solar radiation in the earth's 
atmosphere are water vapor in the troposphere and ozone in the stra- 
tosphere. Water vapor absorbs primarily in the near-infrared region, 
0.7 Mm < A < 4 pm. Ozone (0 3 ) is effective in the ultraviolet (A < 
0.35 Mm) and in the visual (0.5 Mm < A < 0.7 Mm). Fig. 4, taken from 
Pettit (1951) is representative of the spectral absorption for clear 
sky conditions. 

The present parameterization in the model is based on the 
UCLA/GLAS GCM scheme and is described by Davies (1982). It includes a 




Figure 4: Spectral energy curve of solar radiation at sea level and extrapolated outside the 
atmosphere, as given by Pettit (1951). The darkened areas represent gaseous absorption in the 
atmosphere. 
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parameterization for the major absorption processes in the stra- 
tosphere, in the troposphere and at the earth's surface. The parame- 
terization is a function of the water vapor distribution, the cloud 
coverage, the zenith angle of the sun, the albedo of the earth's sur- 
face and the ozone distribution. In this scheme ozone absorption and 
water vapor absorption are assumed to be in above mentioned seperable 
spectral regions. Multiple scattering is taken into account whenever 
it is significant. Following is the brief description of these 
absorption processes. 

a) Water Vapor Absorption 

According to Lacis and Hansen (1974) absorption by water vapor is 
the major source of solar radiative heating in the atmosphere. 
However, parameterizing water vapor absorption is difficult for three 
reasons. These reasons are: 

i. The absorption coefficient of water vapor is highly frequency- 
dependent and accurate monochromatic values are not available for the 
entire spectrum. 

ii. Significant scattering and absorption can occur in the same part 
of the atmosphere . 

iii. The absorption coefficient has a significant dependence on 
pressure. 

For clear skies the effect of' scattering in the spectral region 
of significant water vapor absorption is negligible. The absorption 
due to water vapor is obtained directly from empirical absorption 
functions. The total fractional absorption as a function of water 
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vapor amount is given by (Yamamoto, 1962) 


( y ) - 


2 

(1+141. 5y)°.635 + 5.925y 


where 


y = centimeters of precipitable water vapor. 


(3.85) 


An approximate correction is made for the pressure (p) and 
temperature (T) dependence of the absorption by using an effective 
water vapor amount in place of y as: 


y eff = y(-^-) n (- \ >* 

where n = pressure scaling taken to be 1. 

Pq = reference pressure (1013 mb) 

To = reference temperature (273 K) 


(3.86) 


The absorption by water in the ith layer of the clear atmosphere 
is given by 

Ajj = HqSq {A wv (Myjj + |) - A wv ( My j ) + Rg [A wv (My j ) - A wv (My j+j ) ] } 

’ (3.87) 

Sq = solar influx at the top of the atmosphere 

Mo = cosine of the solar zenith angle (0 q) 

Rg = ground albedo 

y s = effective water vapor amount traversed by the direct solar 
beam in reaching the top of the !th layer 

= I f 1 Q (£__) (To)* d P (3.88) 

g J Po T 

0 

q = specific humidity 

M = magnification factor (Rodgers, 1967) accounting for the 
slant path and refraction. 
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My** 


= 35/ (3 224|i 0 + J)% 

= effective water vapor amount traversed by the reflected 
radiation in reaching the Ith layer from below 


■ g 


i=> ~ g 

= M f q (£_) (T 0 j^ dp + _5_ f q(E_) (T 0 ) % dp 
g J P 0 T" 3g J P 0 T 

0 Pi+1 


(3.89) 


where p* , p* + j and Pg are the pressure at the top of the ith layer, at 
the bottom of the ith layer and at the ground. 

For cloudy sky multiple scattering is the principal factor 
affecting the absorption; in this case the parameterization is based 
on a discrete probability distribution for the absorption coefficient 
derived from measured absorptivities. 


The absorptivity is expressed by: 

00 

A wv (y) = 1 - / p(k)e _k Vdk 


(3.90) 


where 

p(k)dk = fraction of the incident flux that is associated with 
the absorption coefficient between k and k + dk. 

A wv is approximated by the following finite sum 


A wv (y) =1-1 P( k i )e kiy 

i=l 


(3.91) 


where only 5 significant absorption coefficients are considered. 
Davies (1982) has calculated values for p{kj) and kj . These values 


are given in Table 2. 
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Table 2 

Discrete probability distribution of water vapor 
absorption coefficients 
From Davies, 1982 


i 

k i 

-*■ 

P(kj.) 

1 

0.005 

0.107 

2 

0.041 

0.104 

3 

0.416 

0.073 

4 

4.752 

0.044 

5 

72.459 

0.025 
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The radiative transfer problem including multiple scattering for 
cloudy sky is solved as discussed by Davies (1982). The optical pro- 
perties of each layer associated with each are specified by the 
total optical thickness £ of the layer and by the single scattering 
albedo S>. 

c 

Cje.i = f i + ki(y 1+1 - yjE ) (3.92) 

°l,i = Ci/Ci,i 0-93) 

c 

where is the optical thickness due to cloud particles in the ith 
layer. 

b) Ozone Absorption 

According to Lacis and Hansen (1974) the major source of stra- 
tospheric heating is the absorption of solar radiation by ozone. The 
absorption due to ozone occurs at the wavelengths where the Rayleigh 
scattering is significant. But, the absorption by ozone can be 
accurately parameterized because the significant absorption takes 
place very high in the atmosphere where there is little scattering. 
Most of the Rayleigh scattering takes place below the ozone layer, so 
the lower atmosphere acts mainly as a reflecting layer. On the other 
hand, the temperature and pressure dependance of the absorption coef- 
ficient are not very large, and can be neglected if the coefficients 
are chosen for a temperature representative of the ozone layer. 

The absorption due to the ozone amount is different in the weak 
visual bands (Chapius band) versus the ultraviolet bands (Hartley and 
Huggins bands). Figure 5 shows the percent of the solar flux absorbed 



Absorption (*/* Solar Flu: 



Figure 5: Percent of the solar flux absorbed as a function of ozone 
amount for two spectral regions. 
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as function of ozone amount for the two spectra] regions. The ultra- 
violet bands are practically saturated at 0.5 cm of ozone, while the 
visual bands remain nearly proportional to the ozone amount. The 
parameterization of absorption due to ozone must accurately portray 
these different bands. The total fractional ozone absorption as a 
function of ozone amount (X) is expressed as follows: 

A oz< X > = A^(X) + A^ S (X) (3.94) 

where : 

A™ (X) = the fraction of incident solar flux that is absorbed in 
oz 

the ultraviolet bands as a function of ozone amount, 

A VlS (X) = the fraction of incident solar flux that is absorbed in 
oz 

the visual bands as a function of ozone amount. 

Lacis and Hansen (1974) have fitted the frequency -integrated absorption 
curves for ozone to the following expressions. 

A ViS (X) = O JO g llM (3.95) 

02 1.0 + 0.42X + 0 . 000323X 2 

A UV (X) = 1 082X + 0.0658X (3 . 96) 

° Z (1-0 + 138 . 6X) 0 • 805 10 + (103. 6X) 3 


The ozone amount traversed by the direct solar beam in reaching 
the 1th layer of the atmosphere is: 

X £ = ujjM (3.97) 


where : 


Ujj = the amount of ozone in cm in a vertical column above the 1th 
layer, 


M = magnification factor accounting for the slant path and 
refraction, as defined earlier. 

M = 35 

(1224 JJt* + 1) 


(3.98) 
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The ozone path traversed by the diffuse radiation illuminating the ith 
layer from below is: 

X| = u t M + 1.9(u t - u £ ) (3.99) 

where : 

= the total ozone amount in a vertical path above the main 
reflecting layer (i.e., the ground for clear skies or the cloud- 
top for cloudy skies). 

The total absorption of shortwave radiation in ith layer due to ozone 
is therefore: 

A i ,oz ~ A 0 z<X| + l> AozCXjg) } + R(|i 0 ){A 0Z (x?) - A vz (Xj+i ) } ] (3-100) 

where : 

R(|i 0 ) = the albedo of the reflecting region, including the effec- 
tive albedo of the lower atmosphere and the ground reflectivity; 
for cloudy skies it is determined by a multi-layer solution. For 
clear skies it is given by the expression of Lacis and Hansen 
(1974) : 

R( 4 0 ) = R a 0x o ) + [1 - R a (u 0 )](l - Ra)R g /U - Ra R g) (3.101) 

where : 

R a (pt 0 ) = effective albedo of the lower atmosphere. For clear 
skies is given by 0.219/(1.0 + 0.816p o ), 

R g = ground reflectivity, 

R a = the spherical albedo of the reflecting region. For clear 
sky it is given by 0.144. 

R a (p 0 ) and R a were determined by least squares fit. 
c) Earth's Surface Absorption 


The shortwave absorption at earth's surface is treated 
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differently for cloudy and clear skies. 

The surface absorption is divided into two parts. One part 
(Agj) is for the wavelength regions where the absorption coefficient 
of water is significant. Second part (A gj2 ) is for the absorption in 
remaining wavelengths. Approximately 35% of the solar flux is con- 
tained in the spectral regions of significant water vapor absorption. 

A g = A g>1 + A gj 2 (3.102) 

where A g = total surface absorption. 

I . Clear Sky : 

Lacis and Hansen (1974) define the surface absorption for clear sky in 
the spectral regions associated with significant water vapor absorp- 
tion as: 

A g j j = p 0 S 0 [0 . 353 — A^y (My^- ) ] ( 1-R g ) (3.103) 

where : 

M 0 S 0 A wv( M yt) = absorbed direct solar flux by water vapor in 
a vertical column above the ground, 

R g = ground albedo. 

The surface absorption in the spectral regions of negligible water 
vapor absorption is defined as: \ 

___ =* 

A g,2 = ^o s 0 [0.647 - R r (M 0 ) - A 0Z (Mu t )](l - R g )/(1 ~ R r R g) (3.104) 
where : 

PoS 0 A oz ( Mu t) = absorbed direct solar flux by ozone in a vertical 
path above the ground. 

R r (fjt 0 ) = the atmospheric albedo due to Rayleigh scattering, and 
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Is given by: R r (p 0 ) = — 0_u£2 

1 + 6.43 M 0 

" jjc 

R r = the spherical albedo of the Rayleigh atmosphere for illumi- 
nation from below. 

= 0.0685 


II . Cloudy Sky : 


For cloudy sky A g j is obtained by multiplying the total 
transmission at the ground for each value of (absorption coef- 
ficient) by the factor (1-R g ) p (kjj and summing over i for the 
discrete distribution of absorption coefficients given in Table 2. 

Ag 2 is computed from the expression given below where Rayleigh 
scattering is neglected. 


Ag,2 " i*o [0.647 A oz (Mu^ ) ] 


=* 

[1 - R a (M 0 )](l ~ R g )/(1 “ R a R g) 

(3.105) 


where u t is the ozone amount in a vertical path above the highest 

_ ~3 )C 

cloud layer, and R a (|i 0 ) and R a refer to the visual cloud albedo. 


This completes the procedure for calculating the short wave 
absorption in the atmosphere and at the earth’s surface which provides 
the heating rates due to short wave absorption. There are several 
limitations of this scheme. For instance, it ignores the effects of 
aerosols and the spectral nature of the surface albedo. 


3.7 Longwave radiation : 

The longwave radiation scheme of the model is also based on the 
UCLA/GLAS GCM scheme and is described by Harshvardhan and Corsetti 
(1984). The basic equations for the longwave upward and downward 
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fluxes are: 

F clr f <P) = / Ay dP [Bp(T s )T y (p,p s ) + /P Bp(T(p' ) ) dT v (P»P ' ) dpt] 

is d P' (3.106) 

Fclr^P) = / Ay dp[ / P By(T(p')) — ^ - ^p i P ' } dp'] (3.107) 

Pf 

This is for a clear sky case, integrated over the spectral range Ai> . 
Here the individual variables are: 

By(T s ) = blackbody flux at the surface temperature T s , 
p s = surface pressure 
T s = surface temperature, 

T(p' ) = air temperature at pressure p', 
p t = pressure at the top of the atmosphere, 

Ty (p . p ' ) = diffuse transmittance between the levels p and p' , 
v = spectral wavenumber. 

The diffuse transmittance in above equations is not only a func- 
tion of temperature and pressure but it varies very rapidly with wave- 
number u, therefore rendering the computation of fluxes difficult and 
computationally inefficient. In practice, models avoid this problem 
by resorting to certain assumptions regarding the distribution of 
spectral lines in the various absorption bands and constructing 
approximate band models. A fast method of computing the transmittance 
functions over the spectral bands of interest has been developed by 
Chou and Peng (1983) and Chou (1984). 

The thermal infrared spectrum of interest is shown in figure 6. 
The molecular line absorption due to the water vapor stretches from 0 
cm -1 to 3000 cm -1 . Where 0-340 cm -1 , 1380-1900 cm -1 are water vapor 
band centers (regions of strong absorption) and 340-1380 cm - *. 
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1900-3000 cm - * are water vapor wings (regions of moderate absorption) . 
Apart from the molecular line absorption there is continuum absorption 
(e-type) due to water vapor. This e-type absorption is also assumed 
to extend from 0-3000 cm -1 . Absorption due to CO2 is from 580 to 800 
cm -: * where 620 to 720 cnT^ is the band center and 580-620, 
720-800 cnT* are the band wings for CO2. O3 absorption is from 980 to 
1100 cm -1 . The spectral ranges 540-800 cm~^ and 980-1100 cm -1 are 
known as 15 (Jim band and 9.6 pm band respectively. Spectral range 
0-3000 cm - * , excluding the 15 pm band and the 9.6 pm band, is termed as 
water vapor band. 


An alternate form of the flux equations is obtained through 
integration by parts of (3.106) and (3.107). 

Fclrt(p) = APB(T(p)) + G(p,p s ,T s ) - G(p,p s ,T(p s ) ) + 

J T(Ps) 3G(P,P' ,T(p' )) dT(p') (3.108) 

T(p) 3T 


F clrl(p) = APB(T(p)) - G(p,p t .T(p t )) + 


f 


T(pt) 3G(p.p , ,T(p’)I f m P n 


(P) 3T 


where : 


B(T, (p) ) = _1_ / B v (T)dv , 

AP AP 

G(p,p\T) = / Av Ty(p,p' )By(T)dP , 

jp ( n n 1 t ^ 3 By ( T ) 

aG (P'P *11 _ r T / D — dp 

3T Jap T v(P,P } 3T 


(3.109) 

(3.110) 

(3.111) 

(3.112) 


and T(p s ) = air temperature at the surface pressure p s . 

In the above r (p.p) is set to 1 by the definition of transmit- 
If the spectral width Ap of the band is sufficiently narrow 


tance . 



B V (T) is replaced by B(T) , a mean value of the planck function. The 
equation for G becomes: 
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G(p,p',T) = B(T) / T y (p,p')du 


(3.113) 


The radiative cooling rate, or the divergence of net flux is the final 
output of the radiation scheme. It is given by: 



The model uses the above equations to calculate the longwave 
radiative cooling/heating rate of the atmosphere. Special con- 
siderations are given for cloudy skies, either overcast or partial 
cloud cover. The longwave surface flux is also calculated. 


To evaluate the G functions for the fluxes given by (3.108, 109) 
the knowledge of the transmittance function is required. The diffuse 
transmittance function associated with molecular line absorption bet- 
ween the pressure levels p and p ' at the wave number v can be defined 
as (Chou, 1984). 

1 

Ty(p.p') = 2 | exp (-dy(p,p' )/ m) p dp (3.115) 

0 

where p = cos0; 6 is the zenith angle 

and d y is the optical thickness given 

P 

by dy(p.p') = | c(p") ky(p",T(p")) dpi (3.116) 

P' 

where c is the absorber concentration and k is the absorption coef- 
ficient which is a function of wavenumber , pressure and temperature. 
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For computational efficiency k is taken to be the absorption coef- 
ficient at a reference pressure p r and temperature T r scaled by a 
factor depending on pressure and temperature. 

=MPr> T r) f (P.T) (3.117) 

therefore equation (3.116) can be written as 

dy(P.P') = k v (p r ,T r ) w(p.p') = d y (w) (3.118) 

P 

where w(p.p’) = f c(p") f(p",T(p")) dp^ (3.119) 

J g 

P 1 

With the above definitions the molecular diffuse transmittance can be 
written as 

T„(p,p') = Ty (w) (3.120) 

The scaling factor f(p,T) for different absorbers in different 

spectral bands are specified by Chou (1984) as: 

f(p.T) = (E_) m exp [r(T-T r )] for water vapor and C0 2 molecular 

p r absorption. (3.121) 

f(p.T) = 1 for Og molecular absorption. (3.122) 

Where parameters m, r, reference temperature (T r ) and reference 
pressure (p r ) are also taken from Chou (1984). For water vapor 
absorption these parameters are given in table 3 and for C0 2 absorp- 
tion they are given in table 4. 

For water vapor, diffuse transmittance associated with e-type 
absorption is given by 

1 

Ty(p,p') = 2 / exp [~<r v (T 0 )u(p,p' )]p dp (3.123) 

0 

where a v is the e-type absorption coefficient and is given by Roberts 
et al. (1976). 


T 0 = 296°K 
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Table 3 

Parameters for the absorption due to 
water vapor molecular lines 
From Chou, 1984 



h 2 o 

band-center 

h 2 o 

band-wing 

15 Jim 
band 

9.6 Jim 
band 

Spectral range 
(cm -1 ) 

0 - 340 

1380 - 1900 

340 - 540 
800 - 980 

1100 - 1380 

1900 - 3000 

540 - 800 

980 - 1100 

P r (mb) 

275 

550 

550 

550 

T r (K) 

225 

256 

256 

256 

HR" 1 ) 

0.005 

0.016 

0.016 

0.016 

m 

1 

1 

1 

1 



Xl H 
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Table 4 

Parameters for the band-center region and the 
band-wing regions for C0 2 (15 Mm band) 


Band wings 


Band-center Narrow Wide 


AM (AM -1 ) 

Pr (mb) 

Pc (mb) 
m 
n 
a 

r (K) 

r(K'l) 


620 - 720 

30 

1 

0.85 

0.56 

3.1 

15.1 

240 

0.0089 


580 - 620 
720 - 760 
300 
1 

0.85 

0.55 

0.08 

0.9 

240 

0.025 


540 - 620 
720 - 800 
300 
1 

0.50 

0.57 

0.04 

0.9 

240 

0.025 
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u is the scaled water vapor amount and is given by: 


P' 

u(P.P' ) = f e(p") exp 


1800 


(- 


T(p" ) 296 


1_) ] q(p") * d 2l (3.124) 

96 J g 


where e is the water vapor pressure and q is water vapor mixing ratio. 

Similar to equation (3.120) for molecular diffuse transmittance 
we can write transmittance for e-type absorption as 

Ty(p,p') = T v (u) (3.125) 

For water vapor, broad transmission function is given by r y (w)T v (u) s 
t u (w,u) . 

With the above expressions of the transmittance functions the function 

G and its derivative 3G/3T required to estimate the IR fluxes can be 

computed. The expression for G given by equation 3.113 is 

G = B(T) / Ty(p,p')dV (3.126) 

Ay 


G can be rewritten as 

G = B(T) / Ty (z)dy or G = B(T) r(z) (3.127) 
Au 

where t(z) = / r y (z)dy (3.128) 

Ay 

Ty(z) is Ty(w) for molecular absorption and r y (u) for e-type absorp- 
tion. In the regions of overlapping absorption due to different 
absorbers or overlapping molecular and e-type absorption Ty(z) is the 
product of the respective transmittance functions. 

Similarly 3G is given by 3G = 3B(T) t(z) (3.129) 

3T 3T 3T 

Given the absorber amounts t(z) for a given band Au can be evaluated 
using line by calculation with the knowledge of the respective 
absorber coefficients. Using such line by line calculations some 
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regression equations can be developed to provide transmittance T as a 
function of absorber amounts. In the model the following regression 
relations for T(z) for the various bands are used: 

a) Water vapor band 

In the band centers the molecular absorption is much stronger 
than e-type absorption, therefore e-type absorption is neglected for 
band centers . For the band wings both types of absorptions are 
included. 

t(w;T) = t(w;250) [1 + oc(w) (T-250) + j3(w)*( T-250) 2 ] (3.130) 

for band centers 

T(w,u;T) = t(w,u;250) [1 + oc(w,u) (T-250) + j3(w, u)* (T-250) 2 ] (3.131) 

for band wings 

where 

<x(w, u) ,oc(w) = regression coefficients, listed in Table 5, 

/5(w,u), 0(w) = regression coefficients, listed in Table 6, 

and 

r(w, u ; 250) , t(w; 250) = standardized transmission functions, listed in 

Table 7. 

The G function for the water vapor band can be evaluated as 

GwPB = Bj^pc ( T ) T ( w ; T ) + Bty/pw (T)t(w,u;T) (3.132) 

where B wpc (T) and B WPW (T) are respectively the average Plank functions 
over the band centers and band wings. 

b) 15 um band 

In the 15 pm band C0 2 , water vapor molecular absorption and water 
vapor e-type absorption are present. The transmittance in this region 


is given by 
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Table 5 The regrewon coefficient, 10* X a(w, u\ for computing the tnmmmson functions 
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Table 6 The regression coefficient, 10* X 0(w, u\ for computing the transmission functions 
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Table- 7 The 

“Planck weighted” transmission function, 10 J X t(w, tr, 250 K), in the water vapor spectral regions. 
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T(w,u,w 1 ) = t(w)t(u)t(w 1 ) (3.133) 

where w, u and Wj are scaled water vapor amount for molecular absorp- 
tion, scaled water vapor amount for e-type absorption and scaled CO 2 
amount respectively. 

t(w) = exp [-6 . 7w(l + 1 6w° ■ 6 ) ] (3.134) 

t(u) = exp (-2741T 83 ) (3.135) 

t(wj ) = exp (-aw 1 /(l + bw x n ) ) (3.136) 

where a, b and n have different values for CO 2 band center and wings. 
These values are listed in table 4. The final value of z(w^) for the 
entire band is taken to be frequency weighted average of the transmit- 
tence t c (w 1 ) in the band center region and t w (wj) in the band wing 
region, i.e. 

t(W]_) = _1_ [t c (w 1 )Ap c + T^WjJAi^] (3.137) 

Al> 

where Av c , Ai> w and Av are width of the band center, band wing and 
entire 15 pm band respectively. 

The function G is given by 

Gl5/im = B 15Mm (T)T(w)T(u)T( Wl ) (3.138) 

where B 15flm (T) is the average Plank function over the 15 pm band. 

9.6 pm band 

In the 9.6 pm band the water molecular absorption is neglected. 
Thus the total absorption is due to ozone and water vapor e-type 
absorption. Therefore the transmittance for this band is given by 

t(u,w 2 ) = t(u)t(w 2 ) (3.139) 

where w 2 is the ozone amount and u is the scaled water vapor 
amount for e-type absorption. 


r(u) = exp (-9 . 79u) 


(3.140) 
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and 


t(w 2 ) = 1 


83 .31 
120.0 


1 - expj-4.398p ( (1 + 4*345.28 w 2 )^ - 1)}] 

— 0787 p — J 


(3.141) 


P2 P 2 

where p = f (0 3 concentration) *dp/ f qdp 


Pi Pi 

q is the mixing ratio of water vapor. 

The corresponding G function for this band is given by 

G 9-6pm = B 9-6(im (T)T(u)T{w 2 ) (3.142) 

where B 9 g ^(T) is the average Plank function over 9.6 pm band. 


Finally the Plank function in equations 3.108 and 3.109 for 

upward and downward longwave radiation flux is given by 

B(T(p)) = B^pc(T(p)) + B^p^(T(p) ) + B 1 5j im (T(p)) + B 9 g| Ira (T(p)) 

(3.143) 


The function G is given by 

G(P.P'.T) = G wpb + Gjsjju,, + G 9 epm (3.144) 

A similar expression for 3G can also be obtained. 

3T 

The upward and downward flux defined earlier are only valid in 
the clear sky. To account for the effect of cloudy sky the equations 
are modified as described below. 

For longwave radiation the clouds are assumed to be non- 
reflective. For cloudy sky the modified equations (3.108) and (3.109) 
for the fluxes are given by: 

F cldy f (P) = AVB(T(p)) + c(p,p s ) [G(p,p s ,T s ) - G(p,p s ,T(p s )] 

T(p s ) 

+ f c(p,p') 3G (P.P’ ,T(p' ))dT(p' ) 

J 3T 

T(p) 


(3.145) 
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F cldy ; (P) = Ai>B(T(p) ) - c(p,p t ) G(p,p t> T(p t )) 

T(p t ) 

+ f C(p.p') 3G (p,p' ,T(p'))dT(p ! ) (3.146) 

J 3T 

T(p) 

where c (p , p ' ) Is the probability for clear line of sight from p to p ' . 
c (p , P 1 ) depends upon the kind of cloud overlap present between p and 
P ' • For the downward flux c (p ,p ' ) between two pressure levels p and 
P 1 ( P>P ' ) designated by vertical indices i and j respectively is given 
by 

c i j = U-Nj) (1-N j+1 ) (1-Ni.j) (3.147) 

where Nj is the fractinal cloudiness at the level j. N.- lies between 
0 and 1. For clear sky Nj is zero. Similarly for the upward flux 
between the pressure levels p and p' (p’>p) designated by the vertical 
indicies i and j c(p,p') is given by 

Cij = (1-Nj ) (l-Nj.j) (1-N i+1 ) (3.148) 

After determining the probabilities for clear line of sight cjj's the 
radiative fluxes are calculated and then longwave radiation cooling 
rates are computed by equation (3.114). Finally both shortwave and 
longwave cooling rates are added to the thermodynamic equation. 

3 . 8 Surface Energy Balance : 

Over the oceans the prescribed sea surface temperature is used as 
the earth's surface temperature. Over land areas the solution of a 
surface energy balance equation is used to obtain the ground tem- 
perature Tg. The surface energy balance is expressed by, 


( s w4 ~ s wt) + " L wt) ~ F.t F Lt “ 0 


(3.149) 
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i.e. . 

net shortwave flux + net longwave flux - sensible heat flux - latent 
heat flux = 0. In this balance we ignore the flux into the ground and 
the heat capacity of the ground is assumed to be zero. If a denotes 
the surface albedo we express the above equation by, 

(1-a) S W j + (L W | - L W | ) - (Fg-f + F^-t ) = 0 (3.150) 
The surface albedo is specified from climatological tabulation of 
Posey and Clapp (1964) and Kondratyev (1972). Figure 7a illustrates 
the field of albedo for the summer season. We furthermore express the 
fluxes by the following expressions, 


^wt " a Tg 4 


F st = p c H c p |v 2 | (e 2 - e g ) 

F Lt = P C q L |V 2 | ( q 2 _ qg ) GW 


(3.151) 

(3.152) 

(3.153) 


The above were discussed in the surface similarity theory. 

The energy balance equation is coupled to the surface similarity 
theory and is solved by the Newtonian iteration scheme. 


3.9 Orography : 

The orography was extracted from a basic Navy data tape. The 
Navy tape contains global orography on a 10 x 10 minute area. This 
data was averaged at the respective resolution. 


Another field of data that enters the regional model is the sea 
surface temperature (SST). Climatological mean in a current 10 day 
averaged SST are used to define this field. Fig. lh illustrates a 







CHAPTER 4 


PREDICTION OF THE LANDFALL OF A TROPICAL STORM 


During May 1979 a major tropical cyclone formed in the southern 
Bay of Bengal and made a landfall over the east coast of southern 
India around the 12 ^h of that month. This storm caused major destruc- 
tion in life and property over Andhra Predesh. Fig. 8, based on the 
Guam summaries, illustrates the observed track of this tropical 
cyclone . 

The formative aspects of this tropical cyclone have been studied 
by Low Nam (1982). This storm formed along the lower tropospheric 
cyclonic shear zone of the Bay of Bengal monsoonal flows, that was 
well prior to the onset of monsoon over India. The formative period 
of this cyclone coincides with the onset of monsoon and heavy rain over 
Burma. An explosive growth of zonal flow accelerations in the mon- 
soonal low level flows and rather marked divergent outflows over the 
southern Bay of Bengal were noted during the formative period. 

The study, presented here, addresses the landfall of this tropi- 
cal cyclone. The horizontal resolution is 1° lat/long grid. The multi- 
level regional model was initialized with the final FGGE IHb data 
sets for May 11 1979, 00 UTC. At this stage the tropical cyclone 

was located at 13° N and 82.5° E over the Bay of Bengal. 
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Figure 8: Observed track of the tropical cyclone based on the Guam summaries. 
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The ECMWF analysis describes a weaker storm than a typical tropical 
cyclone. The inner rain area of a tropical cyclone is usually con- 
fined to an area within the 150 km radius from the storm's center. 
The strongest winds are nearly always located within this radius from 
a tropical cyclone. At a horizontal resolution of 50 or 100 km an 
adequate description of the inner rain area is not possible since this 
region has to be described by only a few grid points. The ECMWF ana- 
lysis incorporates high resolution cloud drift winds from the geosta- 
tionary satellite (GOES Indian ocean). Those winds were almost the 
only observations over the open oceans in the initial data set. The 
horizontal resolution of these cloud winds were of the order of one 
wind every 200 km 2 . Thus the initial data and the final ECMWF analy- 
sis were not capable of describing the structure of inner rain area 
adequately. In spite of the data deficiencies and coarse horizontal 
resolution the information content of the storm environment 
was somewhat reasonable. The model thus seems to describe a tropical 
storm rather than a tropical cyclone. It is the landfall of such a 
weaker disturbance that is possible to model in the present context. 

4 • * Observational aspects of the landfall 

Departing from a normal representation of flow fields on pressure 
surfaces we shall illustrate the horizontal motion field on the so 
called sigma surfaces. We simply want to look at results for one 
example of a forecast which are made on sigma surfaces without con- 
verting the results to the pressure surfaces, some of which lie below 
the earth's surface. The sigma to pressure conversion does smooth the 
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results to some extent and it is of considerable interest to see these 
uncontaminated results. Here we shall show the circulations at o = 
0.85, o = 0.5 and a = 0.25 surface, these closely correspond to the 
850, 500 and 250 mb surfaces. 

Fig. 9 (a, b, c, d, e, f and g) illustrates the horizontal motion 
field on the <7 = 0.85 surface based on the observed wind analysis. The 
salient feature in these charts is the storm whose westward passage 
starting on May 11 1979 00 UTC through May 14 1979 00 UTC is 
illustrated here in steps of every 12 hours. The storm’s maximum, 
based on the ECMWF analysis remains close to 20 ms- 1 for the first 84 
hours, thereafter a weakening in its intensity was noted as the storm 
passed over south-central India. The velocity maximum remains to the 
east and south-east of the storm center during its westward passage. 
A large scale counterclockwise circulation to the south of this storm 
denotes the southern-equatorial trough whose center was located close 
to 90E and 10S during this entire 4 day period. Southwesterly flows 
extended from the central Bay of Bengal towards Burma throughout this 
period and the Burmese monsoon had become fairly active. The low 
level flows over the Arabian sea were mostly from the North- 
Northwesterly direction, the dry Arabian air is a characteristic 
feature of the pre-onset (over India). 

Of interest to the proposed prediction experiment is the westward 
motion and the landfall of the tropical cyclone which is reasonably 
described by these observations. The circulations based on the ECMWF 
analysis over the middle troposphere are described in the sequence of 
charts, Fig. 10 (a, b, c, d, e, f and g). These illustrations show 




Figure 9 a, b, c, and d: Horizontal motion field on the a = 0.85 surface based on the observed wind 
analysis starting May 11 1979 00 UTC to May 12 1979 1200 UTC. 
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the horizontal motion field on the a = 0.5 surface for the period May 
11 1979 00 UTC through May 14 1979 00 UTC in steps of 12 hours. 

A major feature at the a = 0.5 surface is the interaction of a so 
called westerly trough with the storm’s upper level circulation. At 
the initial time this trough passes through 35N and 60E. It shows a 
steady eastward propagation and is located near 70E and 20N by the 
14 th D f May. Subsequently this westerly trough axis merges with the 
upper circulation of the tropical cyclone. During this later stage 
considerable dry air enters the storm circulation as it weakens. The 
wind speed at the a = 0.5 surface is generally less that 20 ms- 1 , 
however the cyclonic circulation is quite well marked. A ridge 
extends northwards over Burma where the warm monsoonal air resides at 
this time. 

The warm region is where the upper tropospheric anticyclonic 
circulation is found, which can be seen from the charts at the a = 0.25 
level after the initial date (see fig. 11a through g) . This is the 
monsoonal upper anticyclone which makes its way to the Tibetan Plateau 
later in the summer season. The westerly upper trough which was noted 
at the a = 0.5 surface has a stronger intensity at this level (<x = 
0.25). Ahead of this upper trough strong winds of the order of 56 
ms - 1 were noted initially. These strong winds move well ahead of the 
trough in the subsequent 4-day period. A closer inspection of the 
flows over the tropical cyclone reveals a pronounced upper level 
diffluent geometry. That feature was present during this entire 
sequence. A part of the flow over the tropical cyclone moves northward 
and eventually eastwards around the upper anticyclone. Another part 
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Figure H a, b, c and d: Horizontal motion field on the a = 0.25 surface based on the observed wind 
analysis starting May 11 1979 00 UTC to May 12 1979 1200 UTC. 
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of the same flow moves westwards and eventually southwards around an 
upper cyclonic circulation located over southern India. These two 
features i.e., the upper anticyclone over Burma (and Indochina) and 
the cyclonic circulation over southern India appear to be closely 
related to the upper circulation of the tropical cyclone. 

It is of interest to see if several of the aforementioned 
features can be predicted by the regional model. 

4 . 2 Prediction of the landfall 

The predicted motion field showed very little error in the first 
72 hours. The predicted flow fields closely paralleled those based on 
observations. We shall first illustrate the predicted fields at the 
a = 0.85 surface (close to 850 mb). A sequence of predicted charts at 
hours 12, 24, 36, 48, 60 and 72 are shown in fig. 12 (a, b, c, d, e 
and f) respectively. Nearly all of the circulations features are 
reasonably placed in the prediction. The landfall of the tropical 
cyclone occurred 36 hours after the initial time along the southwest 
coast of India. The model predictions at hour 36 places the center of 
the storm's circulation at almost the same position on the east coast 
of India as was noted from the observed fields in fig. 9. The sub- 
sequent positions of the storm after landfall at hours 48, 60 and 72 
show a remarkable agreement with the observed track. At hour 72 the 
observed center at a = 0.85 was located at roughly 79E and 18N while 
the predicted position was roughly at 79E and 19N. When the storm was 
located over land (at hour 72) the strongest winds were located 
offshore. The observed (i.e., based on ECMWF analysis) wind maxima 
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a 


b 


c 


Figure 12 a, b and c: Predicted horizontal motion field on the a = 0.85 
surface starting May 11 1979 1200 UTC to May 12 1979 1200 UTC. 
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had an intensity around 15 ms- 1 , the model results show stronger winds 
of the order of 24 ms- 1 . The model winds over the ocean were in 
general stronger. This discrepancy in the intensity of winds may be 
due to a systematic error of the model. Alternatively, it is possible 
that the cloud drift winds were not assigned a correct level, thus the 
analysis may have a bias. The current numerical model does not 
include the vertical flux of momentum by shallow or deep cumulus 
clouds, which is another factor that is worthy of consideration. 

At the middle troposphere (o = 0.5) the storm (based on observed 
analysis) has a somewhat stronger intensity and is located somewhat 
westwards of the position at a = 0.85. A westward tilt with height 
has been noted for such disturbance after landfall by Douglas 1987. 
The tilt is related to the increase of westerly wind with height in 
such hydrostatic systems. The model prediction at a = 0.5 level are 
illustrated in fig. 13 (a, b, c, d, e and f) at intervals of every 12 
hours. Again we note a remarkable success in the prediction of cir- 
culations through the 3 day period. The position of the storm at a = 
0.5 is nearly identical in the observed and the predicted panels at 
hour 72. The predicted winds are however weaker than the observed; 
around the storm. Thus it appears that the model prediction of the 
wind speeds are somewhat overestimated at a = 0.85 and are underesti- 
mated at ff = 0.5. To the east of the storm the speed of the southerly 
flow, for both the observed and the predicted, are around 20 ms- 1 . 
The major discrepancy is to the west of the storm where the model’s 
underestimation is clearly apparent at the a - 0.5 level. This regi 
over land does show the presence of many shallow and towering cumuli. 
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It Is possible that the addition of cumulus momentum flux in the model 
might rectify this bias. 

At the a = 0.25 level, fig. 14 (a, b, c, d, e, and f), the model 
handles the eastward movement of the westerly trough remarkably well. 
The model also handles the strong winds ahead of the westerly trough 
quite well; the model's maximum winds over this region lie between 40 
to 50 ms- 1 during the first two days, thereafter the model exhibits a 
further strengthening of the winds reaching almost 66 ms- 1 by 72 
hours. Over this region near 35N and 75E, the observed winds at a = 

0.25 level were in excess of 50 ms- 1 . The predicted winds at the 

jet stream level seem to be somewhat too high - although this region 
of the western Tibetan Plateau does not have adequate upper air 
observations. The upper anticyclone over Burma is very well repre- 
sented by the predictions. The strong westerly winds to its north 
around 310E and 37N compare very closely - the observed maximum wind 
is around 70 ms- 1 while the predicted value is around 69 ms- 1 . 

Overall, the predicted flow fields at all of these levels show a 

marked agreement with observations. 

4.3 Rainfall prediction : 

The total precipitation in the model comes from the parameteriza- 
tion of cumulus convection and from the disposition of super- 
saturation. The latter is here identified as the non-convective 
contribution. 


A recent atlas of daily rainfall rates, Krishnamurti et al. 
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(1983c), describes a procedure for the estimation of rainfall rates 
from a mix of satellite and raingauge observations. The FGGE lie data 
collection includes a very large archive of raingauge observations. 
Well over 3,000 observations/day provide daily rainfall totals over a 
MONEX domain ( 30E to 150E, 30S to 40N) . In view of these practical 
limitations, TIROS N infrared radiance data was used. This is a polar 
orbiting sun synchronous satellite which makes roughly 2 local passes 
a day with an equator crossing around 3AM and 3PM local time. The 
raingauge data are available on 24 hourly intervals at 00 UTC 
each day. The method of analysis consisted in first performing a 
statistical regression among the satellite radiances, its time rate of 
change and the raingauge rainfall totals over 24 hourly periods. For 
this purpose all of the satellite observations were color ated onto the 
raingauge sites via a space-time interpolation. The time synchroniza- 
tion required the use of cubic splines, and thus the time rate of 
change of radiances for this resolution represents an integrated 
representation of the growth (or decay) of large cloudy areas. 
Although such a use of data sets from polar orbiting satellites has 
not been made previously, it appears quite promising for NWP because 
of the space-time scales. The next step in this analysis is the 
construction of a first guess rainfall field for objective analysis 
from the aforementioned regression coefficients. This is done using 
the daily grid point values of the radiances and their time rate of 
change. The next step is an objective analysis of the raingauge rain- 
fall field. This makes full use of all available raingauge data sets 


from the FGGE lie data. 
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A comparison of the 24 hour precipitation amounts based on the 
satellite and raingauge with those obtained from the two versions of 
the model are presented in figs. 15, 16 and 17 (a, b, c) . In each of 
these figures panel (a) shows the ’observed' 24 hourly total and the 
other two panels show the predicted estimates. In all of these charts 
the interval of analysis is 10 mm/ day . The respective totals in figs. 
15, 16 and 17 (a, b, c) are for the hours 0-24, 24-48 and 48-72. The 
panels (b) in these illustrations are based on a terrain tabulation 
that was obtained on a 2.5° latitude/longitude grid and interpolated 
on to the model resolution of 0.9375° latitude/longitude. The panels 
(c) illustrate the results from using terrain tabulations from a high 
resolution U.S. Navy data set on the resolution of the model. In 
general the model-based rain fall rates during the landfall of the 
tropical cyclone are weaker and spread out over a larger horizontal 
area. This discrepancy is perhaps related to the resolution. 
Although the arrival of heavy rain is reasonably predicted the posi- 
tion errors are still somewhat too large. In the first 24 hours the 
'observed' coastal rainfall over Tamil Nadu was quite large compared 
to the predicted estimates. The model places a center of heavy rain 
over the open ocean. Although the circulation forecasts seemed fairly 
accurate during this period. A number of factors can contribute to 
this discrepancy. We believe that the initial humidity analysis over 
the open oceans is based entirely on the analysis-forecast cycle of 
ECMWF . Use of satellite radiance data to improve the humidity analy- 
sis following Krishnamurti et al. (1984) and Cadet (1983) can provide 
some improvement. The model prediction overestimates the strong 
southerly flows to the east of the tropical cyclone. That may also 
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have been an important factor for the discrepancy in the rainfall over 
the ocean. The cyclonic shear region is shifted eastward, although 
the circulation center is reasonably placed. The stronger convection 
on the cyclonic shear side of the predicted flows is a consistent com- 
putational result. This region also has stronger mass and moisture 
convergence. Better resolution and initialization procedures, such as 
dynamic relaxation (Nudging), may be needed to overcome some of these 
problems. The forecasts of the rainfall rates for hours 24-48 and 
48-72 are quite reasonable, they show the landfall of a hurricane with 
heavy rain occurring over southern India, however as before the 
quality of forecast (district by district) still requires much impro- 
vement. The use of a finer resolution terrain (panels (c) of figs. 
15, 16 and 17) appears to shift the rainfall maximum somewhat south- 
wards initially when the storm center is over the open ocean. As the 
storm approached the Bay of Bengal Coast the interaction of the lower 
level circulations with the eastern Ghats results in an enhancement of 
rainfall total along the northern coasts of the peninsular. 



CHAPTER 5 


PARAMETER I ZATON OF GROUND WETNESS 


The interaction between the atmosphere, the land surface and the 
biosphere plays an important role in determining the weather through 
radiation, heat, moisture and momentum exchanges. For more than a 
decade, evidence has been accumulating from climate and forest model 
experiments that the atmosphere is sensitive to variations in pro- 
cesses at the land surface. The experiments of Charney (1975) and 
Charney et aJ[. (1977) demonstrated considerable sensitivity to albedo 
and the parameterization of evaporation respectively. One essential 
factor in the surface evaporation estimation is the parameterization of 
soil moisture. Many soil moisture sensitivity experiments (Walker and 
Rountree, 1977; Shukla and Mintz, 1982; Yeh et al . , 1984; Sud and 
Fennessy, 1984; Sud and Smith, 1985b; Kitoh et al . , 1988; Sud and 
Molod, 1988) have confirmed the importance of its parameterization in 
climate models. The general conclusion of these experiments is that a 
reduced soil moisture causes a reduction in evaporation which in turn 
decreases rainfall. Accurate initial conditions were shown to be 
important in the context of predictions on timescales from a few hours 
or days (Walker and Rountree, 1977; Rountree and Bolton, 1983) up to 
several months (Carson and Sangster, 1981). It is evidently important 
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to represent land surface processes realistically; this requires an 
adequate parameterization scheme, including specification of suitable 
data sets of vegetation, soil type and initial soil moisture. 

At the earth’s surface available (net radiative) energy is par- 
titioned between three fluxes: 1. Surface soil heat flux, 2. 

Evaporative or latent heat flux and 3. Sensible heat flux. The 
balance of these fluxes with the net radiative fluxes determines the 
changes in surface temperature. To estimate the evaporative flux the 
knowledge of available soil moisture is necessary. Soil moisture also 
has a large impact on the thermal properties of land surface such as 
heat capacity and thermal conductivity. 

Soil moisture/Ground wetness is a difficult parameter to estimate 
as it depends on the evaporation, rainfall, snowmelt, infiltration, 
surface runoff as well as vegetation. For a layer of ground between 
depths z and z + 6z, neglecting horizontal sub-surface transfer of 
water the generalized equation for soil moisture change can be written 
as 

am s (z) = _ mm + N (z) ( 5 -n 

3t 3z 

where m s (z) = soil moisture at depth z 
M(z) = moisture flux at depth z 
and N(z) = source/sink term. 

At the surface, moisture flux is given as 

M(0) = -E(0) + (P + M s - Y(0) ) (5.2) 

where E = evapotranspiration 
P = precipitation rate 
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M s = snowmelt 
and Y = surface runoff. 

Ground wetness parameter can be defined as fractional soil 
moisture at the surface: 

GW = m s (Q) (5.3) 

mS max 

where m s max is the maximum value (field capacity) of m s (0). 

To parameterize the above many simplifications are made in the 
present day models.' One of the widely used methods in the climate 
models was first introduced by Manabe (1969a) resulting from Soviet 
observational studies by Romanova (1954). Romanova found that for 
plains and forest regions most of the moisture change takes place in 
the top one meter of soil, which usually encompasses the root zone of 
moist vegetation. Manabe’ s scheme uses one layer water budget 
equation which includes a simple parameterization for evaporation, 
precipitation, surface runoff and snowmelt. Field capacity is taken 
constant for all land points. The ground wetness sensitivity experi- 
ments mentioned earlier and others reviewed by Mintz (1984), use rela- 
tively simple parameterization schemes and lack geographical 
variations of most surface characteristics. More complex schemes have 
been developed in recent years which include the geographic variations 
of albedo, roughness length and root depth (Hansen jet al . , 1983; 

Dickinson, 1984; Sellers et al. , 1986). 

All leading numerical regional models use very simple parameteri- 
zation for the ground wetness. NCAR/Penn State limited area model, 
JMA Tokyo regional spectral model and ECMWF limited area model use 
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climatological specified values for ground wetness. French weather 
service limited area model uses one layer predictive equation for soil 
moisture. In the present formulation of FSU regional model we use one 
of the following empirical relations for ground wetness. 


GWM1 


( a max a ) 
^“max " “min) 


(1 - RH) 


(5.4) 


or 


GWM2 = 0.85 [1.0 - exp{-200.0(0.25 - a) 2 }] (5.5) 

where a is surface albedo and RH is relative humidity for the nearest 
layer to the ground. Climatological specified values of ground wet- 
ness or given by some empirical relations like equations (5.4) and 
'5. 5) do not provide the change in soil moisture due to rainfall and 
other moisture fluxes at the surface. If a predictive equation is 
used for GW (Ground Wetness) as in French weather service model, one 
faces the problem of getting reliable initial values for the ground 
wetness besides the difficulties of parameterizing the various fluxes 
involved. 


In the present chapter we propose a new scheme for ground wetness 
parameterization. The procedure can be used to obtain the initial 
values for the ground wetness to initialize the model. Further, it 
provides an expression for ground wetness which depends on the rain- 
fall and other large scale surface parameters. The following section 
describes the proposed scheme. In the last section 5.3 some results 
of the new scheme are presented and the results are compared with the 


old scheme. 
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5.1 The Parameterization Scheme : 

The scheme is based on the moisture budget analysis described by 
Haiyan He jet <il . (1987), First, the surface evaporation is estimated 
from vertically integrated apparent moisture sink (defined later in 
this section) and the rainfall rates. Using this evaporation rate 
surface energy balance is solved for the surface temperature. 
Thereafter surface evaporation and temperature provide ground wetness 
using similarity theory. 
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w ^ ere ^ _ latent heat of vaporization 

Ps = surface pressure 

PT = pressure at the top of the column taken as 125 mb 
Q 2 = apparent moisture sink 

= ~ L (la + V»Vq + u iS) (5.10) 

3t 3p 

R = Precipitation rate per unit area at the surface 
Fqs = Surface evaporation per unit area 

= Fq(p s ) 

F Q(Pt) = vertical moisture flux at the top of the atmosphere 
which is assumed to be zero. 

We define LFq S = Surface latent heat flux 

Q 

s F L 

Equation (5.9) can be rewritten to provide surface latent heat 
flux as: 

Q Ps 

F L = LR - 1 / Q 2 dp (5.11) 

g Pt 

Equation (5.11) is used to calculate F L . For computing Q 2 
observed data for horizontal wind components (u,v) and specific 
humidity FGGE level II B data is used. The vertical p-velocity w is 
obtained from the horizontal divergence by integrating the continuity 
equation 

3u + 3v + 3u = 0 (5.12) 

3x 3y 3p 

with the surface boundary condition on u at p = p s expressed by: 

w = w s = -g p s ( u s M + v s IS) (5.13) 

a~cos$ 3X a dip 

This boundary condition implies that the surface flow is parallel 
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to the surface 


terrain. 


H = H(X,(j>) surface terrain, where X is the longitude and <j> is the 
latitude . 


u s 1 v s = surface zonal and meridional component of V 
Assuming that the motion is approximately adiabatic near the tro- 
popause, an additional condition is imposed on the vertical velocity oi 
at p T (s: 125 mb) from the thermodynamic equation i.e.: 


— ^ 96 + V*V8^ 

« = «! = 9t (5.14) 

86 

3p 

where 6 is the potential temperature. Using these constraints the 
horizontal divergence is corrected as: 


D corrected _ D original + 


w T - - / Ps D °rigina] dp 

Pp 


(5.15) 


After obtaining F^, we make use of the surface energy balance and 
similarity theory to compute ground wetness parameter (GW). To 
improve the accuracy of the solution we adapt a soil-slab model for 
surface energy balance. Following Blackadar (1979) let us consider a 
thin surface soil layer at the surface temperature T g Q in thermal con- 
tact with the deeper layer, which is assumed to be acting as a heat 
reservoir at constant temperature T m : 


C g 8T g Q = S w l(l-a) + (L w 4 - <xT g Q) 


F s t " F L t 


- C g K ffl (T g Q - T m ) (5.16) 

The first right hand side term of the above equation is the net short 
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wave radiation flux into the ground. The second term is the net 
longwave radiation flux into the ground. Third and fourth terms are 
the upward sensible heat and latent heat fluxes. The last term provi- 
des the flux due to conduction between the upper thin soil layer and 
the lower substrata. 

a = surface albedo 

Cg and K m are given by the following expressions 

Cg = 0.95 (Ac s /2<»))^ (5.17) 

K m = 1.180) (5.18) 

where 

c s = Heat capacity of the soil per unit volume 
= PgC sg 

c sg is the specific heat capacity of the soil and p g is the soil 
density. 

o) = angular velocity of the earth's rotation 

X = thermal conductivity of the soil. 

Heat capacity of the soil c s and thermal conductivity of the soil 
X are taken as functions of the moisture present in the soil. 

For the present formulation, we have assumed c s and X to be 
linear functions of the soil moisture/ground wetness parameter GW as: 

C s = (1.42 + GW*1 .68)*10 6 Jm _3 k _1 (5.19) 

X = (0.25 + GW*1 . 33) Wnr’k -1 (5.20) 

F s (surface sensible heat flux) in equation (5.16) is calculated 
using the following expression 

F s t = pC H C p |V 2 | (T 2 - T g Q) (5.21) 

where T 2 is the temperature at a chosen level within the constant flux 
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layer at the surface. V 2 is the horizontal wind vector at the same 
level. Cp is the specific heat. Cg is calculated using similarity 
theory discussed in section 3.4. The expression for Cg in stable and 
unstable surface layers are given by equations (3.53) and (3.62). From 
section 3.4 it is evident that for calculating Cg knowledge of surface 
temperature is required. Therefore equation (5.16) is solved iterati- 
vely along with the similarity theory. S K i and L w i are calculated 

using the radiation parameterization described in chapter 3. For Fg 

Q 

in equation (5.16) F^ obtained by the expression (5.11) is used. 
Since c s and X depend on the ground wetness, climatological values for 
GW are used as a first guess. For the subsequent iterations GW is 
updated using the procedure given below. 

Once the surface temperature (TgQ) is obtained surface saturation 
specific humidity q gs is computed using Teten's formula. Using this 
surface temperature stability dependent Cg and Cq are calculated and 
ground wetness parameter GW is estimated using the following formula: 


gwQ = _1_ 

Qgs 


Q2 " 


(5.22) 


PCqL |v 2 | J 

This expression for GW is obtained by writing latent heat flux Fl 
similar to the expression for F s given by equation (5.21). 


F L = PC Q L|V 2 | (q 2 - GW*qg S ) (5.23) 

Q 

The above equation is solved for GW and F L is replaced by F L given by 
equation (5.11). This generates a large data base for the ground wet- 
ness from the vertically integrated moisture budget. This GW is next 
expressed as a function of various large scale variables. A selective 
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second order regression scheme is used to provide GW as a function of 

average 24 hour accumulated rainfall rate, surface temperature, 

earth's surface albedo and terrain as 
„ t .New _ 

GWT - a 0 + ajR + a 2 T g + a 3 oc' + a 4 H + a 11 R® + a 22 T g z 

+ a 33 oc'2 + a 44 H 2 + a 12 RT g + a 13 R<x' + a 14 Rlf 

+ a 23 TgOc + a 24 TgH + a 34 oc H (5,24) 

The regression coefficients are given in Table 8. In the above 
equation R is the 24 hour rainfall divided by 24. T g is the surface 
temperature and <x and H are scaled albedo and scaled terrain height 
given by: 


K '~ (“max °0 (5.25) 

Kax - ^min) 

and H = (H max - H) (5.26) 

H min) 

where o^max' ^max are the maximum values of albedo and terrain. 
Similarly oc min and H m -[ n are the minimum values of albedo and terrain 
fields. As it is evident from Table 8 GW has a strong dependence on 
rainfall. The dependence on terrain field is mainly to account for 
the runoff due to slope of the ground. Ground wetness is closely 
related to the surface temperature also. 

To carry out the above scheme FGGE data for 11 May 1979 to 16, 
May 1979 and 20 June 1979 to 28 June 1979 for horizontal winds, rela- 
tive humidity, temperature geopotential height at all standard 
pressure levels is used. Daily rainfall rates for the above mentioned 
period are extracted from Summer Monex experiment data. The domain 
chosen is 2.8°-38.4°N and 45.9°-113.4°E. The above mentioned 15 days 
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described above are presented in figures 19-23 in comparison with 
ground wetness (GWM1 , GWM2 ) estimates from expressions (5.4) and 

(5.5) . Figures 24 (a, b, c, d and e) show the average 24 hour rain- 
fall (mm/day) for the five subperiods. The observed measures of rain- 
fall, shown here, were derived from a mix of raingauge and satellite 
based estimates. Comparison of GW maps with the rainfall for the 
corresponding period reveals clearly that GW provided by the new 
scheme is in good agreement with the rainfall distribution, whereas 
GWM1 and GWM2 do not represent the soil moisture well. The ground 
wetness for the five subperiods, clearly illustrates the major defi- 
ciency of the specified ground wetness (shown in the bottom panel of 
each diagram). The ground wetness distribution from GWM2 is rather 
flat and exceeds 80 % over most of India and Southeast Asia. The top 
panel on each of the figures 25-29 shows the latent heat flux F L Q over 
the land calculated from moisture budget analysis using expression 

(5.11) for 5 subperiods. The bottom two panels in these figures are 

S 1 s 2 

the maps of latent heat flux (F^ , F^ ) calculated using (5.4) and 

(5.5) expressions for GW, the similarity theory and surface energy 

balance as presently used in the model, discussed in Chapter 3. 

s new 

The second panels on these figures show F^ which is also calculated 

by similarity dependent formula (5.23) but it uses ground wetness 

given by GwQ and slab model for surface energy balance as discussed in 

s new 

section 5.2. It is clearly seen from the figures that F L is in 
close agreement with F L Q whereas the surface latent heat flux distri- 
butions given by F^ .F^ are not very close to FlQ. The two older ver- 
sions of ground wetness fail to represent the surface flux of latent 
heat. Similarly the next 30-34 figures show the corresponding 
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Figure 19 a, b and c: Average GW parameter maps for the subperiod 
1 (May 11 1979 00 UTC - May 13 1979 12 UTC). It ranges 
from 0 for dry land to 1 for saturated land, (a) calcu- 
lated from moisture budget analysis, (b) and (c) are 
calculated respectively by the expressions given by 
(5.4) and (5.5). 
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Figure 20 a, b and c : Average GW parameter maps for the subperiod 
2 (May 14 1979 00 UTC - May 16 1979 12 UTC) . It ranges 
from 0 for dry land to 1 for saturated land, (a) calcu- 
lated from moisture budget analysis, (b) and (c) are 
calculated respectively by the expressions given by 
(5.4) and (5.5). 
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Figure 21 a, b and c: Average GW parameter maps for the subperiod 
3 (June 20 1979 00 UTC - June 22 1979 12 UTC) . It 
ranges from 0 for dry land to 1 for saturated land. (a) 
calculated from moisture budget analysis, (b) and (c) 
are calculated respectively by the expressions given by 
(5.4) and (5.5). 
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Figure 22 a, b and c: Average GW parameter maps for the subperiod 
4 (June 23 1979 00 UTC - June 25 1979 12 UTC) . It 
ranges from 0 for dry land to 1 for saturated land. (a) 
calculated from moisture budget analysis, (b) and (c) 
are calculated respectively by the expressions given by 
(5.4) and (5.5). 







GW OBTRINED FROM Q2 


H5. 9E 

60. 0E 

70. 0E 

80. 0E 

30 . 0E 

o 

o 

o 

m 

1 13. HE 

38.HN 

1 j , 

m rrr 

v y /■/ i 



- 


H5.9E 60. OE 70. OE 80. OE 90. OE 100. OE 113. HE 

Figure 23 a, b and c: Average GW parameter maps for the subperiod 
5 (June 26 1979 00 UTC - June 28 1979 12 UTC). It 
ranges from O for dry land to 1 for saturated land. (a) 
calculated from moisture budget analysis, (b) and (c) 
are calculated respectively by the expressions given by 
(5.4) and (5.5) . 
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Figure 24 a, b and c: Three-day average observed rainfall rates 
(mm/day) derived from a mix of raingauge and satellite based estimates 
for the period (a) May 11 1979 00 UTC - May 14 1979 00 UTC, (b) May 14 
1979 00 UTC - May 17 1979 00 UTC, (c) June 20 1979 00 UTC - June 23 1979 
00 UTC. 
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Figure 24 d and e: Three-day average observed rainfall rates (mm/day) 
derived from a mix of raingauge and satellite based estimates for the 
period (d) June 23 1979 OO UTC - June 26 1979 00 UTC, (e) June 26 1979 
00 UTC - June 29 1979 00 UTC. 
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Figure 25 a, b, e and d: Average surface latent heat flux (watts/m* 2 ) for the subperiod 1 (May 11 1979 
00 UTC - May 13 1979 12 UTC) calculated from (a) moisture budget analysis, (b) similarity theory and 
GW from moisture budget analysis, (e) similarity theory and GW from expression (5.4), (d) similarity 
theory and GW from expression (5.5). 
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Figure 26 a, b, c and d: Average surface latent heat flux (watts/m 2 ) for the subperiod 2 (May 14 1979 
00 UTC - May 16 1979 12 UTC) calculated from (a) moisture budget analysis, (b) similarity theory and 
GW from moisture budget analysis, (c) similarity theory and GW from expression (5.4), (d) similarity 
theory and GW from expression (5.5). 






137 



Figure 27 a, b, c and d: Average surface latent heat flux (watts/m 2 ) for the subperiod 3 (June 20 
1979 00 UTC - June 22 1979 12 UTC ) calculated from (a) moisture budget analysis, (b) similarity theory 
and GW from moisture budget analysis, (c) similarity theory and GW from expression (5.4), (d) simi- 
larity theory and GW from expression (5.5). 





138 



Figure 28 a, b, c and d: Average surface latent heat flux (watts/m 2 ) for the subperiod 4 (June 23 
1979 00 UTC - May 25 1979 12 UTC) calculated from (a) moisture budget analysis, (b) similarity theory 
and GW from moisture budget analysis, (c) similarity theory and GW from expression (5.4), (d) simi- 
larity theory and GW from expression (5.5). 
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Figure 30 a, b, c and d: Average surface temperature (°K) for the subperiod 1 (May 11 1979 00 UTC - 
May 13 1979 12 UTC). (a) Calculated from moisture budget analysis, (b) calculated using surface 
latent heat flux shown in figure 25 (b), (c) calculated using surface latent heat flux shown in figure 
25 (c), (d) calculated using surface latent heat flux shown in figure 25 (d) . 
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Figure 33 a, b, c and d: Average surface temperature (°K) for the subperiod 4 (June 23 1979 00 UTC - 
June 25 1979 12 UTC). (a) Calculated from moisture budget analysis, (b) calculated using surface 
latent heat flux shown in figure 28 (b), (c) calculated using surface latent heat flux shown in figure 
28 (c), (d) calculated using surface latent heat flux shown in figure 28 (d) . 
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Figure 34 a, b, c and d: Average surface temperature (°K) for the subperiod 5 (June 26 1979 00 UTC - 
June 28 1979 12 UTC). (a) Calculated from moisture budget analysis, (b) calculated using surface 
latent heat flux shown in figure 29 (b), (e) calculated using surface latent heat flux shown in figure 
29 (c), (d) calculated using surface latent heat flux shown in figure 29 (d) . ’ 
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surface temperatures calculated from each of the schemes for 5 sub- 
periods . 

Next, we shall divide the domain of interest into four 

s new S 1 s 2 

subdomains; the averages of latent heat flux (F^Q, F L , f l > F L ) 
over the land are computed for each of these subdomains. These avera- 
ges are compared in figures 35 and 36 for the five subperiods. The 
subdomains are defined as: 


I 

12. 2 C 

- 38 . 4°N, 

45.9° 

- 66 . 6°E 

II 

o 

00 

CM 

- 21 . 6°N, 

66.6° 

- 91 . 9°E 

III 

21.6° 

- 38.4°N, 

66.6° 

- 91 . 9°E 

IV 

12.2° 

- 38 . 4°N, 

co 

CO 

o 

- 113.4° 


The figures 35-36 further confirm that for all the subdomains 
average flux calculated from the new parameterize tion of ground wet- 
ness is very close to the surface flux F L Q determined by large scale 
Q 2 and rainfall rate. The fluxes using GWM3 and GWM2 are not in good 
agreement with F L Q. Finally the figures 37 and 38 present the ground 
wetness (GW^ ew ) given by expression (5.24). This is the final form of 
the ground wetness as function of rainfall, ground temperature, sur- 
face albedo and terrain height. GW New has a correlation coefficient 
of 0.7 with GW Q. 

A 72 hour forecast of the tropical cyclone discussed in Chapter 4 
is made to see the impact of improved parameterization of ground wet- 
ness. The model is initialized for May 11 1979, 00 UTC and a three 
day forecast is made. The figures 39 and 40 show the time evolution 
of the ground wetness as predicted in this experiment. Ground wetness 
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Figure 37 a, b and c': Average GW calculated from the GW expression 
given by (5.24). This is based on the regression analysis of GW estima- 
tes from moisture budget analysis with the parameters past 24 hour rain- 
fall, surface temperature, albedo and terrain. (a) Subperiod 1 (May 11 
1979 00 UTC - May 13 1979 12 UTC), (b) Subperiod 2 (May 14 1979 00 UTC - 
May 16 1979 12 UTC), (c) Subperiod 3 (June 20 1979 00 UTC - June 22 1979 
12 UTC). 
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Figure 38 a and b: Average GW calculated from the GW expression given 
by (5.24). This is based on the regression analysis of GW estimates 
from moisture budget analysis with the parameters past 24 hour rainfall, 
surface temperature, albedo and terrain. (a) Subperiod 4 (June 23 1979 
00 UTC - June 25 1979 12 UTC), (b) Subperiod 5 (June 26 1979 00 UTC - 
June 28 1979 12 UTC). 




Figure 39 a, b, c and d: Maps showing GW at 12 hours interval a? predicted by the model using new 
GW parameterization for (a) May 11 1979 00 UTC, (b) May 11 1979 12 UTC, (c) May 12 1979 00 UTC, and 
(d) May 12 1979 12 UTC. 
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maps are 12 hours apart. By referring to the observed rainfall shown 
earlier in figures 15a, 16a and 17a it can be inferred that soil 
moisture is very well predicted by the new scheme. Next figure 41 
shows the predicted 3-day average surface temperature and latent heat 
flux. A comparison of figures 18-a, 18-c and 41-b shows that average 
latent heat flux using dynamic GW new is enhanced compared to the pre- 
dicted average using specified ground wetness as GWM2 given by 
(5.5). Regions of high latent heat fluxes as seen in figure 18-a are 
better predicted with the new dynamic ground wetness compared to 18-c. 
Figure 42 shows the time series of predicted rainfall and ground wet- 
ness at two grid points located at 31.9 °N, 71.3 °E and 23.4 °N, 79.7 
°E. In figures 43 and 44 time series of surface temperature and 
latent heat flux using old and new parameterization of GW are pre- 
sented at the same grid points. The comparison shows that the diurnal 
changes in surface temperature and Fl are predicted better with the 
new parameterization scheme for GW. Fl with new scheme is enhanced 
for both the points. The results on the diurnal change of the ground 
temperature do not exhibit a sharp step function like behavior when 
the ground wetness is parameterized by the proposed regression 
approach. Two prominent features in these results are i) The ampli- 
tude of the latent heat fluxes from the regression based ground wet- 
ness parameterization are much larger in comparison to the estimate 
from the specific ground wetness: ii) The amplitude of the diurnal 
change is very well preserved cycle after cycle for the regression 
based ground wetness parameterization. The precipitation predictions 
using the new parameterization scheme are presented in the next 


chapter . 
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a 



Figure 41: (a) Average latent heat flux (watts/m 2 ) for May 11 1979 00 UTC 

- May 13 1979 12 UTC predicted by the model using new GW parameteriza- 
tion. Shaded areas show the high values of latent heat flux. (b) 
* Average (May 11 1979 00 UTC - May 13 1979 12 UTC) surface temperature 
(°K) as predicted by the model in the above mentioned experiment. 
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the above mentioned grid points. 




















CHAPTER 6 


SENSITIVITY OF RAINFALL PREDICTION TO GROUND 
WETNESS AND RESOLUTION 


In this Chapter we shall compare the results of forecasts on the 
sensitivity to the ground wetness parameterization, discussed in the 
previous Chapter, for two different horizontal resolutions. 

The specific example we explore here is the same case study that 
was discussed in Chapter 4 on the landfall of a tropical storm over 
southern India. 

Two satellite photographs shown in figs. 45a and 45b illustrate 
the banded structure of the cloud cover of the tropical cyclone just 
prior to and after the landfall. These illustrations were obtained 
from the DMSP satellite. They describe the cloud cover at 10:00 am 
local time over India for May 12 and May 13, 1979, respectively. The 
cloud cover during these two days shows the typical structure of the 
tropical cyclone spiral rainbands. The dense cloud cover of the eye 
wall is also clearly portrayed by these illustrations. An open eye is 
also evident on the satellite picture of May 12th when the storm was 


over the ocean. 



40 N 
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Figure 45a: Satellite photograph showing the cloud cover at 10:00 am local time over India for May 
12, 1979, 
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Figure 45b: Satellite photograph showing the cloud cover at 10:00 am local time over India for May 
13, 1979, 


160 


It Is of considerable interest to see the extent to which such 
features and details are revealed in the predicted rainfall fields. 
We shall next examine the 'observed' raingauge based rainfall records 
for these periods. 

In figures 46a and 46b we show 24 hourly rainfall totals based on 
the FGGE data collections for the respective periods May 11 00 UTC to 
May 12 00 UTC and May 12 00 UTC to May 13 00 UTC. The observed rain- 
fall amounts are derived from roughly 3000 raingauge sites over the 
domain illustrated here. The rainfall amounts were averaged over 1 
degree latitude/longitude squares in order to smooth out local 
variability and to facilitate the illustration of this vast volume of 
data. A careful examination of these observed fields, shown in 
figures 46a and 46b clearly show a banded geometry in the rainfall 
distribution. The highest rainfall amounts are of the order of 100 to 
160 mm/day in the near coastal regions of Madras and Andhra states. 

In the control experiment the ground wetness was simply expressed 
as a function of the prescribed surface albedo. According to this 
formulation, the ground wetness has a strong inverse relationship to 
the surface albedo. In the parameterized ground wetness experiment 
the statistical relation, described in Chapter 5, was used. The 
remaining aspects of the model for these two sets of experiments were 
kept identically the same. Furthermore the data sets and the initial 
states were also identical . In a two to three day prediction 
experiment we do not expect a large impact from the parameterization 
of land surface processes. That is more of a climate issue; in cli- 
mate modelling significant impact of surface wetness parameters have 
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Figure 46: Accumulated 24 hour observed rainfall derived from raingauge 
data. Raingauge data has been averaged over 1 degree lat/long squares, 
(a) For the period May 11 1979 00 UTC to May 12 1979 00 UTC (b) For the 
period May 12 1979 00 UTC to May 13 1979 00 UTC. 
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been demonstrated by many Investigators such as: Yamakazi (1989), 
Rind (1984) and several others mentioned in the previous chapter. 
Basically these studies show a major improvement in the climatological 
rainfall patterns with a time dependent modelling of the ground wet- 
ness. A reduced soil moisture causes a reduction in evaporation which 
in turn is shown to reduce the rainfall amounts in these studies. A 
motivation for carrying out the proposed present comparison study is 
to see if the parameterization scheme is stable and if it shows any 
measurable changes in the precipitation estimates within the con- 
sidered time frame. 

The other aspect of experimentation covered in this section is on 
the combined effects of resolution and ground wetness parameterization 
on the prediction. The increase of resolution does have a direct 
effect on the increase of vertical motion and rainfall rates. The 
increased rainfall provides an enhancement of the ground wetness para- 
meterization and a positive feedback via latent heat fluxes towards 
further increase of rainfall can be expected. 

The spiral rainbands in a hurricane have been studied for a 
number of years. In the context of numerical weather prediction 
models the simulation of the convection along the spiral rainbands and 
along the eye wall of a hurricane is recognized as a problem of major 
interest . 

Three dimensional hurricane forecasting (or simulation) studies 
began with the efforts of Anthes et al . (1971) and Miller et al ■ 
(1972). In recent years Yamasaki (1986 and 1989) has reported on some 
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6 • 1 Rainfall prediction at 1° latitude/longitude resolution : 

The predicted rainfall rates for the first 24 hours and second 24 
hours for 1° latitude/longitude resolution are respectively presented 
in figures (47 a,b) and (48 a,b). In both the figures, panel (a) is 
predicted rainfall from the forecast using simple ground wetness para- 
meterization (based on surface albedo) and panel (b) shows predicted 
rainfall from the forecast using improved ground wetness parameteriza- 
tion (as discussed in chapter 4). By comparing the figures (47 a,b 
and 48 a,b) with the corresponding observed rainfall rates presented 
in figures (46 a,b) the sensitivity of the rainfall prediction to the 
ground wetness can be observed. It is clearly seen that rainfall 
rates over the land are enhanced in the new ground wetness experiment. 
For the first 24 hour rainfall, the rainfall over Bangladesh, 
South-East coast of India and Burma-Thailand coast is predicted better 
in the forecast using new ground wetness scheme. Similarly for the 
second 24 hours rainfall the position and intensity of rainfall maxima 
over land as shown by the shaded areas in these figures are better 
predicted by the new ground wetness forecast. 

6 . 2 Rainfall prediction at 0.5° latitude/longitude resolution : 

Figures (49 a,b) and (50 a,b) show respectively the first 24 
hours and second 24 hours predicted rainfall for a forecast at hori- 
zontal resolution 0.5° latitude/longitude. Panel (a) in these 
diagrams shows the results for the simple representation of ground 
wetness while panel (b) shows the results for the improved ground wet- 
ness parameterization. At this horizontal resolution of the model 



165 



4-5 . 9t0 , 7 SQ.GE 2 , 2 7C * 3£ 80. OE 30. 0E 100. OE 113. 4E 


2H HO URS TOTAL _RRINFnLL_ENDING 00GHXJ 2 19791 NEW GW RUN } INT*10 



Figure 47: Accumulated 24 hour predicted rainfall for the period May 11 
1979 00 UTC to May 12 1979 00 UTC. Model resolution is .9375° lat/long. 
(a) From a control experiment using old GW parameterization (b) From a 
sensitivity experiment using new GW parameterization. 
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Figure 48: Accumulated 24 hour predicted rainfall for the period May 12 
1979 00 UTC to May 13 1979 00 UTC. Model resolution is .9375° 
lat/long. (a) From a control experiment using old GW parameterization 
(b) From a sensitivity experiment using new GW parameterization. 
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Figure 49: Accumulated 24 hour predicted rainfall for the period May 11 
1979 00 UTC to May 12 1979 00 UTC. Model resolution Is .46875° 
lat/long. Rainbands are highlighted by solid lines, (a) From a control 
experiment using old GW parameterization (b) from a sensitivity experi- 
ment using new GW parameterization. 
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Figure 50: Accumulated 24 hour predicted rainfall for the period May 12 
1979 00 UTC to May 13 1979 00 UTC. Model resolution is .46875° 

lat/long. Rainbands are highlighted by solid lines, (a) From a control 
experiment using old GW parameterization (b) From a sensitivity experi 
ment using new GW parameterization. 
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both the experiments produce rainbands. However, there is an 
interesting and important positive impact from the improved ground 
wetness parameterization. The 0-24 and. 24-48 hours rainfall rates 
shown in these figures bear a very strong resemblance to the banded 
structures of the satellite photographs shown in figures 45 (a,b). 
The location of the central areas of the spiralling predicted rainfall 
bands and the satellite cloud cover both converge to the same general 
storm center. The error in the location of the storm center is quite 
small at both the 24 and 48 hour panels in particular for the improved 
ground wetness parameterization experiment. A number of details of 
banded structure over south eastern Arabian sea, peninsular India, Bay 
of Bengal and Malay-Indonesia regions are very well represented. As 
stated earlier there are, however, some major differences in the pre- 
dictions with the two representations of the ground wetness parame- 
terization. 

The differences between the rainfall estimates for the two ground 
wetness parameterization experiments are presented in figures 51 (a) 
and (b) . Panel (a) shows the difference for first 24 hours rainfall 
and panel (b) shows the difference for second 24 hours rainfall. The 
shaded areas in these figures show the regions where the rainfall 
rates predicted by the new ground wetness experiment are higher com- 
pared to the control model run (using simple parameterization for GW) . 
The banded structure of the rainfall is enhanced in the new ground 
wetness experiment. The rainfall over land is in general higher with 
the new GW scheme. In particular some major improvements in the 
second day rainfall prediction over Southern India with the new GW 
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Figure 51: The difference in 24 hour predicted rainfall between model 
forecasts using old and new GW parameterization for .46875" lat/long 
resolution for the period (a) May 11 1979 00 UTC to May 12 1979 OO UTC 
and (b) May 12 1979 00 UTC to May 13 1979 00 UTC. The shaded regions 
indicate rainfall with new GW > rainfall with old GW. 
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scheme are shown in figures 52 (a) and (b). These figures are the 
same as figures 50 (a) and (b) but the contouring interval in these 
diagrams is 25 mm/day. On the South-East coast of India around 80° E 
and 13°N very heavy rainfall was observed as seen in figure 46 (b) for 
the period May 12 1979 00 UTC to May 13, 1979 00 UTC. The arrival of 
this heavy rainfall is very well predicted by the forecast using new 
GW scheme as seen in figure 52 (b). The position and magnitude of 
rainfall center over land is in very good agreement with the observed 
rainfall. In the control experiment there is a rainfall maximum over 
the ocean located at 82°E, 13°N and the rainfall over land around 
80°E, 13°N is of the order of 50 mm/day, which is unde predicted. 
Similarly the rainfall over South-Central India along the 78°E longi- 
tudinal line is very well predicted by the forecast using new GW 
scheme. The maximum rainfall along this line over Southern India is 
of the order of 75 mm/day from this forecast which is the same as 
raingauge observed rainfall. Whereas in the control run maximum 
rainfall along this line is of the order of 50 mm/day. Overall the 
precipitation pattern shows a closer agreement to the observed 
raingauge estimates for the prediction with the new ground wetness 
parameterization. 

Finally, it is observed from these results that the impact of 
ground wetness parameterization is more on the second day rainfall 
forecasts. This confirms the fact that a realistic parameterization 
of ground wetness is very important for rainfall forecasts on larger 


time scales. 
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Figure 52: Accumulated 24 hour predicted rainfall for the period May 12 
1979 00 UTC to May 13 1979 00 UTC. Model resolution is .46875° 

lat/long. (a) From a control experiment using old GW parameterization 
(bl From a sensitivity experiment using new GW parameterization. This, 
figure is same as figure 50 a, b but the contouring interval is 25 

mm/day , 




CHAPTER 7 


CONCLUSION 


The most important results of this study are (1) the numerical 
prediction of the detailed rainbands of tropical cyclone as it made 
its landfall over India and (2) the improvement of rainfall rates pre- 
diction over the land using the new ground wetness parameterization 
scheme. The mechanism for the formation of rainbands are not 
addressed in this thesis. The high resolution model advects the 
cyclonic storm to nearly its correct position over southern India and 
the bands are nearly correctly located as compared to rather complex 
bands from a high resolution satellite photograph. These results per- 
tain to roughly 48 hours after the initial state. The horizontal 
resolution of the model was roughly 50 km. At a resolution of 
approximately 100 km, these banded features were not adequately 
resolved. The impact of improved ground wetness was observed to be 
more on the second day rainfall prediction. The overall conclusion of 
the forecast experiment presented in chapter 4 is that the circulation 
forecasts made by the model are of quite high accuracy. However the 
details of the precipitation forecasts are reasonably good on the 
synoptic scale but require further improvements on the meso-scale. 
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The success of these results have a bearing on many aspects of the 
model development that are addressed in this thesis. 

• The semi-Lagrangi an advective scheme is shown to be a power- 
ful method for the handling of the nonlinear advective terms. In the 
context of the shallow water equation this scheme has been shown to 
conserve the parcel and domain invariants very well to within 1% for 
periods of the order of 5 days. 

• The computation time was saved by a factor of 4 to 5 from the 
use of a semi-implicit time differencing scheme. In comparison, the 
use of explicit time differencing, takes considerably more computer 
resources . 

• The use of Arakawa c-grid with a staggering in the horizontal 
has enabled us to simplify the lateral boundary conditions for the 
pressure functions. 

• The solution of a three dimensional Helmholtz equation was 
simplified from the use of a similarity transformation. This has 
enabled us to transform this problem into a series of two dimensional 
Helmholtz equations, these were solved from the use of a fast tri- 
diagonal matrix solver. 

• The use of dynamic normal model initialization was found to 
adequately describe the initial divergence and precipitation fields. 
This scheme developed by Sugi (1981) seems to converge quite rapidly 
at the highest horizontal resolution (*= 50 km) presented here. 
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The use of steep mountains, called 'envelope orography' was 
incorporated into the prediction model. This has enabled the incor- 
poration of steep mountains such as the Himalayas, the western Ghats 
and other orographic features around the Indian subcontinent. 

• The modified Kuo's scheme, as used in the present study, was 
shown to handle precipitation rates that well exceed 100 mm/day. 

• Although no sensitivity studies on the radiative transfer 
formulations are presented here, the model does include a state-of- 
the-art band model. This band model evaluates the transfer of short 
and longwave irradiances in the presence of modelled clouds, diurnal 
change and surface energy balance. 

• An area where considerable research was carried out is the 
area of surface flux calculations of latent heat. Here a parameteri- 
zation procedure for the ground wetness parameter was first developed. 
This required diagnostic calculations of moisture budget from 
atmospheric observations. Vertical integral of the moisture budget 
provided a data base for the surface flux of latent heat. These 
fluxes were incorporated into surface energy balance. Thereafter 
using surface similarity theory a data base for the ground wetness 
parameter was generated. The data for the ground wetness was next 
statistically regressed against a number of prediction model parame- 
ters and a parameterization of the ground wetness parameter was 
derived. This parameterization was next used in prediction experi- 
ments and compared with more conventional parameterization of the 
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ground wetness. The results of these experiments provided the best 
results at high resolution and from the use of the parameterized 
ground wetness parameter. 

• The use of the parameterized ground wetness parameter was 
shown to provide a reasonable diurnal cycle of the ground temperature 
and surface fluxes. The later were in close agreement with values 
obtained from observed moisture budget estimates. 

This study has several limitations. Perhaps the most obvious one 
is the lack of use of predicted boundary conditions. The specified 
boundary conditions used here are so called 'perfect values' based on 
observation. This seems to provide the best results. The forecasts 
were seen to degrade from the use of time invariant or global model 
forecast based boundary values. Further work is needed in this area. 
The model is quite sensitive to the use of various physical parame- 
terizations such as cumulus convection and the planetary boundary 
layer. The ground wetness parameterization developed in this work 
does not make use of the vegetation cover over land. Further detailed 
sensitivity studies are required to provide a better understanding in 
these areas. In order to address the predictability over the monsoon 
region a larger number of experiments are needed. 
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